script_rwtBPDN_iterative

 BPDN rwt update (initialize with standard homotopy)

 Solve the following basis pursuit denoising (BPDN) problem
 min_x  \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2

 Initialize with unweighted BPDN
 and dynamically update the weights w_i

 Written by: Salman Asif, Georgia Tech
 Email: sasif@ece.gatech.edu
 Created: March 16, 2011

0001 % BPDN rwt update (initialize with standard homotopy)
0002 %
0003 % Solve the following basis pursuit denoising (BPDN) problem
0004 % min_x  \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2
0005 %
0006 % Initialize with unweighted BPDN
0007 % and dynamically update the weights w_i
0008 %
0009 % Written by: Salman Asif, Georgia Tech
0010 % Email: sasif@ece.gatech.edu
0011 % Created: March 16, 2011
0012 
0013 function out = script_rwtBPDN_iterative(in);
0014 A = in.A; y = in.y; x = in.x;
0015 gamma_orig = find(x);
0016 [M N] = size(A);
0017 tau = in.tau; max_rwt = in.max_rwt;
0018 x_init = in.x_init;
0019 
0020 %% parameter selection
0021 % max_rwt = 3;
0022 % lambda = 1e-2;
0023 % tau = min(lambda*max(abs(A'*y)),sigma*sqrt(log(N)));
0024 
0025 rwt_mode = in.rwt_mode; 
0026 
0027 delx_mode = in.delx_mode;
0028 homotopy_update = 'v1';
0029 
0030 maxiter = 2*N;
0031 iter_ALL = [];
0032 err_ALL = [];
0033 time_ALL = [];
0034 supp_diff = [];
0035 
0036 if norm(x_init) == 0
0037     %% Standard BPDN
0038     in = [];
0039     in.tau = tau;
0040     in.maxiter = maxiter;
0041     in.x_orig = x;
0042     in.record = 1;
0043     in.delx_mode = delx_mode;
0044     tic;
0045     out = BPDN_homotopy_function(A, y, in); %BPDN
0046     % time_ALL = [time_ALL out.time];
0047     time_ALL = [time_ALL toc];
0048     
0049     x_init = out.x_out;
0050     gamma_old = out.gamma;
0051     iter_old = out.iter;
0052     iter_ALL = [iter_ALL iter_old];
0053     err_ALL = [err_ALL norm(x-x_init)/norm(x)];
0054     
0055     supp_diff = [supp_diff length(setxor(gamma_old,gamma_orig))];
0056 end
0057 W_new = ones(N,1)*tau; % homotopy solved with same fixed weights.
0058 
0059 %% Iterative reweighting
0060 xh_mod = x_init;
0061 gamma_xh = gamma_old;
0062 for rwt_itr = 1:max_rwt
0063     gamma_old = gamma_xh;
0064     x_old = xh_mod;
0065     W_old = W_new;
0066     
0067     [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_old,M);
0068     W_new = tau/alpha./(beta*abs(x_old)+epsilon);
0069     
0070     % W_new = tau/epsilon*ones(N,1);
0071     % W_new(gamma_old) = min([tau*ones(length(gamma_old),1) tau./(beta*abs(x_old(gamma_old)))],[],2);
0072     
0073     switch homotopy_update
0074         case 'v1'
0075             pk_old = A'*(A*x_old-y);
0076             pk_old(gamma_old) = sign(pk_old(gamma_old)).*W_old(gamma_old);
0077             in = [];
0078             in.x_old = x_old;
0079             in.gamma = gamma_old;
0080             in.pk_old = pk_old;
0081             in.W_old = W_old;
0082             in.W_new = W_new;
0083             dW = W_new-W_old;
0084             in.maxiter = maxiter;
0085             
0086             % delx = -AtAgx\(-dW(gamma_old).*sign(pk_old(gamma_old)));
0087             % in.delx = delx;
0088             in.delx_mode = delx_mode;
0089             switch in.delx_mode
0090                 case 'mil';
0091                     % in.delx_mode = 'mil';
0092                     % The following gram matrix and its inverse can be used from the
0093                     % previous homotopy. Too lazy to include that right now...
0094                     % wt BPDN homotopy update
0095                     AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0096                     iAtAgx = inv(AtAgx);
0097                     in.AtA = AtAgx;
0098                     in.iAtA = iAtAgx;
0099                 case {'qr','chol'};
0100                     % in.delx_mode = 'qr';
0101                     [Q R] = qr(A(:,gamma_old),0);
0102                     in.Q = Q; in.R = R;
0103                 case 'qrM'
0104                     % in.delx_mode = 'qrM';
0105                     [Q0 R0] = qr(A(:,gamma_old));
0106                     in.Q0 = Q0; in.R0 = R0;
0107             end
0108             tic
0109             out = wtBPDN_Update_function_v1(A, y, in);
0110             % time_update = out.time;
0111             time_update = toc;
0112             xh_mod = out.x_out;
0113             gamma_xh = out.gamma;
0114             iter_update = out.iter;
0115             
0116         case 'v2'
0117             
0118             % Homotopy update v2 (weighting embedded in the matrix A)
0119             AW_old = A*diag(tau./W_old);
0120             u0_hat = x_old.*(W_old/tau);
0121             ds = AW_old'*(AW_old*u0_hat-y);
0122             AW = A*diag(tau./W_new);
0123             yhat = AW*u0_hat;
0124             pk_old = ds; pk_old(gamma_old) = sign(pk_old(gamma_old)).*tau;
0125             in = [];
0126             in.x_old = u0_hat;
0127             in.gamma = gamma_old;
0128             in.pk_old = pk_old;
0129             in.tau = tau;
0130             in.maxiter = maxiter;
0131             in.yhat = yhat;
0132             
0133             % ds = pk_old; d = AW'*(yhat-y);
0134             % delx = AtAgx\(ds(gamma_old)-d(gamma_old));
0135             % in.delx = delx;
0136             in.delx_mode = delx_mode;
0137             switch in.delx_mode
0138                 case 'mil';
0139                     % in.delx_mode = 'mil';
0140                     AtAgx = AW(:,gamma_old)'*AW(:,gamma_old);
0141                     iAtAgx = inv(AtAgx);                    
0142                     in.AtA = AtAgx;
0143                     in.iAtA = iAtAgx;
0144                 case {'qr','chol'};
0145                     % in.delx_mode = 'qr';
0146                     [Q R] = qr(AW(:,gamma_old),0);
0147                     in.Q = Q; in.R = R;
0148                 case 'qrM'
0149                     % in.delx_mode = 'qrM';
0150                     [Q0 R0] = qr(AW(:,gamma_old));
0151                     in.Q0 = Q0; in.R0 = R0;
0152             end
0153             tic
0154             out = wtBPDN_Update_function_v2(AW, y, in);
0155             time_update = out.time;
0156             time_update = toc;
0157             xh_mod = out.x_out.*(tau./W_new);
0158             gamma_xh = out.gamma;
0159             iter_update = out.iter;
0160     end
0161     iter_ALL = [iter_ALL iter_update];
0162     err_ALL = [err_ALL norm(x-xh_mod)/norm(x)];
0163     time_ALL = [time_ALL time_update];
0164     supp_diff = [supp_diff length(setxor(gamma_xh,gamma_orig))];    
0165     
0166     %     % compare with solution from scratch
0167     %     in = [];
0168     %     in.tau = tau;
0169     %     in.maxiter = maxiter;
0170     %     in.x_orig = x;
0171     %     in.record = 1;
0172     %     AW = A*diag(tau./W_new);
0173     %     out_new = BPDN_homotopy_function(AW, y, in); %BPDN
0174     %     xh_mod = out_new.x_out.*(tau./W_new);
0175     %     gamma_new = out_new.gamma;
0176     %     new_iter = out_new.iter;
0177     %     new_time = out_new.time;
0178     
0179     %     figure(1);
0180     %     subplot(221); plot([x xh xh_v2 xh_mod]);
0181     %     subplot(222); plot([xh-xh_v2 xh-xh_mod]);
0182     %     subplot(223); plot(xh,x,'.');
0183     %     subplot(224); plot([W_new W_old]);
0184     %     err_rec = norm(xh-x)/norm(x);
0185     %     err_old = norm(x_init-x)/norm(x);
0186     % disp(sprintf(['RWT iter. = %d, wtBPDN-LASSO homotopy iter. = %d, Update = %d, Update-v2=%d, Rec. error = %3.4g, improve = %3.4g'], rwt_itr, new_iter, iter_update, update_iter_v2, err_rec , (err_old-err_rec)/err_old*100));
0187 end
0188 
0189 out.x_out = xh_mod;
0190 out.x_init = x_init;
0191 out.iter = iter_ALL;
0192 out.gamma = gamma_xh;
0193 out.W_new = W_new;
0194 out.err = err_ALL;
0195 out.time = time_ALL;
0196 out.supp_diff = supp_diff;