l1homotopy_v1

 l1homotopy_v1.m

 A general program that solves homotopy for
   a weighted LASSO/BPDN problem with or without a warm-start vector
   
   Some examples for dynamic updating include
   sequential measurements
   time-varying signal
   iterative reweighting
   measurement replacement
   dictionary learning
   Kalman-type filtering
   Streaming signal recovery 

   or any other problem that can be written in the following form:

 minimize_x  \Sum [1-epsilon]w_old + epsilon w_new]|x_i| + 1/2*||Ax-(1-epsilon)yh-epsilon y||_2^2 + (1-epsilon)u'x,

 where solution is updating by changing epsilon from 0 toward 1.
   yh = A*xh_old, 
   u = -w.*sign(xh_old) (where w can be the old or the new weights)

 Optimality conditions:

 A'(Ax-yh)+u + epsilon(A'(yh-y)-u) = -(w_old + epsilon(w_new-w_old))z    on Gamma
 |A'(Ax-yh)+u + epsilon(A'(yh-y)-u)| < w_old + epsilon(w_new-w_old)    off Gamma

   
 Inputs:
  A -- M x N measurement matrix
  y -- measurements

  opts - input structure

   xh_old  -- old signal estimate
   gamma   -- support of xh_old

   (LHS terms)
   pk_old  -- A'(A*xh_old-yh)+u
   Atdyu   -- A'(yh-y)-u

   (RHS terms)
   W_old  -- w_old
   W_new  -- w_new

   AtAgx and iAtAgx (i.e., A(:,gamma)'*A(:,gamma) and its inverse)
   or QR/Cholesky factors

   delx_mode -- mode for rank-1 update ('mil', 'chol', or 'qr')
   nonneg  -- add nonneg constraint on the solution? (default = 0)
   maxiter -- maximum number of homotopy iterations
   Te      -- maximum support size allowed
   record  -- record iteration history
   x_orig  -- origianl signal for error history
   debias  -- debias the solution at the end
   early_terminate -- terminate early if the support is identified
                   (useful only in high SNR settings)
   verbose -- print output after every verbose steps
   plots   -- plot the solution at every iteration after verbose steps

 Outputs:
 out - output structure
   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
   error_table - error table with iteration record
   iAtA on Gamma, or QR/Cholesky factors

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

-------------------------------------------+
 Copyright (c) 2012.  M. Salman Asif
-------------------------------------------+

0001 % l1homotopy_v1.m
0002 %
0003 % A general program that solves homotopy for
0004 %   a weighted LASSO/BPDN problem with or without a warm-start vector
0005 %
0006 %   Some examples for dynamic updating include
0007 %   sequential measurements
0008 %   time-varying signal
0009 %   iterative reweighting
0010 %   measurement replacement
0011 %   dictionary learning
0012 %   Kalman-type filtering
0013 %   Streaming signal recovery
0014 %
0015 %   or any other problem that can be written in the following form:
0016 %
0017 % minimize_x  \Sum [1-epsilon]w_old + epsilon w_new]|x_i| + 1/2*||Ax-(1-epsilon)yh-epsilon y||_2^2 + (1-epsilon)u'x,
0018 %
0019 % where solution is updating by changing epsilon from 0 toward 1.
0020 %   yh = A*xh_old,
0021 %   u = -w.*sign(xh_old) (where w can be the old or the new weights)
0022 %
0023 % Optimality conditions:
0024 %
0025 % A'(Ax-yh)+u + epsilon(A'(yh-y)-u) = -(w_old + epsilon(w_new-w_old))z    on Gamma
0026 % |A'(Ax-yh)+u + epsilon(A'(yh-y)-u)| < w_old + epsilon(w_new-w_old)    off Gamma
0027 %
0028 %
0029 % Inputs:
0030 %  A -- M x N measurement matrix
0031 %  y -- measurements
0032 %
0033 %  opts - input structure
0034 %
0035 %   xh_old  -- old signal estimate
0036 %   gamma   -- support of xh_old
0037 %
0038 %   (LHS terms)
0039 %   pk_old  -- A'(A*xh_old-yh)+u
0040 %   Atdyu   -- A'(yh-y)-u
0041 %
0042 %   (RHS terms)
0043 %   W_old  -- w_old
0044 %   W_new  -- w_new
0045 %
0046 %   AtAgx and iAtAgx (i.e., A(:,gamma)'*A(:,gamma) and its inverse)
0047 %   or QR/Cholesky factors
0048 %
0049 %   delx_mode -- mode for rank-1 update ('mil', 'chol', or 'qr')
0050 %   nonneg  -- add nonneg constraint on the solution? (default = 0)
0051 %   maxiter -- maximum number of homotopy iterations
0052 %   Te      -- maximum support size allowed
0053 %   record  -- record iteration history
0054 %   x_orig  -- origianl signal for error history
0055 %   debias  -- debias the solution at the end
0056 %   early_terminate -- terminate early if the support is identified
0057 %                   (useful only in high SNR settings)
0058 %   verbose -- print output after every verbose steps
0059 %   plots   -- plot the solution at every iteration after verbose steps
0060 %
0061 % Outputs:
0062 % out - output structure
0063 %   x_out - output for BPDN
0064 %   gamma - support of the solution
0065 %   iter - number of homotopy iterations taken by the solver
0066 %   time - time taken by the solver
0067 %   error_table - error table with iteration record
0068 %   iAtA on Gamma, or QR/Cholesky factors
0069 %
0070 % Written by: Salman Asif, Georgia Tech
0071 % Email: sasif@gatech.edu
0072 %
0073 %-------------------------------------------+
0074 % Copyright (c) 2012.  M. Salman Asif
0075 %-------------------------------------------+
0076 
0077 function out = l1homotopy_v1(A, y, opts)
0078 
0079 N = size(A,2);
0080 M = size(A,1);
0081 
0082 % Use prior, related information if provided
0083 % old solution
0084 if isfield(opts,'xh_old'); xh_old = opts.xh_old; else xh_old = zeros(N,1); end
0085 % old solution constraints
0086 if isfield(opts,'pk_old'); pk_old = opts.pk_old; else pk_old = -A'*y; end
0087 % old solution support
0088 if isfield(opts,'gamma'); gamma_xh = opts.gamma; else gamma_xh = []; end
0089 if isempty(gamma_xh); [tau_old gamma_xh] = max(abs(pk_old)); end
0090 % dummy variable (..)
0091 if isfield(opts,'Atdyu'); Atdyu = opts.Atdyu; else Atdyu = 0; end
0092 % new regularization parameters
0093 if isfield(opts,'W_new'); W_new = opts.W_new; else W_new = opts.tau; end
0094 % old regularization parameters
0095 if isfield(opts,'W_old'); W_old = opts.W_old; else W_old = tau_old; end
0096 % no initial solution
0097 if norm(xh_old) == 0; pk_old = -A'*y; Atdyu = 0; [tau_old gamma_xh] = max(abs(pk_old)); W_old = tau_old; end
0098 % input is a zero vector
0099 if norm(y) == 0; 
0100     out = opts;
0101     out.x_out = zeros(N,1);
0102     out.gamma = []; % find(abs(xk_1)>0);
0103     out.iter = 0;
0104     out.time = 0;
0105     disp('input is a zero vector');
0106     return;
0107 end
0108 % output is a zero vector
0109 if nnz(abs(A'*y) < W_new) == N
0110     out = opts;
0111     out.x_out = zeros(N,1);
0112     out.gamma = []; % find(abs(xk_1)>0);
0113     out.iter = 0;
0114     out.time = 0;
0115     disp('output is a zero vector');
0116     return;
0117 end
0118 
0119 % Make vectors out of scalar regularization parameters..
0120 W_old = ones(N,1).*W_old; W_new = ones(N,1).*W_new;
0121 
0122 % maximum iterations
0123 if isfield(opts,'maxiter'); maxiter = opts.maxiter; else maxiter = 2*N; end
0124 % non-negativity constraint
0125 if isfield(opts,'nonneg'); nonneg = opts.nonneg; else nonneg = 0; end
0126 % maximum support size
0127 if isfield(opts,'Te'); Te = opts.Te; else Te = inf; end
0128 % record error/history
0129 if isfield(opts,'record'); err_record = opts.record; else err_record = 0; end
0130 if err_record; err_fun = opts.err_fun; end 
0131 % debiasing step at the end (solve LS on the support)
0132 if isfield(opts,'debias'); debias = opts.debias; else debias = 0; end
0133 % early terminate if residual of restricted LS falls below certain
0134 % threshold
0135 if isfield(opts,'early_terminate'); early_terminate = opts.early_terminate; else early_terminate = 0; end
0136 % print output
0137 if isfield(opts,'verbose'); verbose = opts.verbose; else verbose = 0; end
0138 % debug plots
0139 if isfield(opts,'plots'); plots = opts.plots; else plots = 0; end
0140 
0141 %% GO
0142 t0 = cputime;
0143 
0144 %% Initial step
0145 epsilon = 0;
0146 Wk = (1-epsilon)*W_old + epsilon*W_new;
0147 dW = W_new-W_old;
0148 
0149 temp_gamma = zeros(N,1);
0150 temp_gamma(gamma_xh) = gamma_xh;
0151 gamma_xc = find([1:N]' ~= temp_gamma);
0152 
0153 z_x = zeros(N,1);
0154 z_x(gamma_xh) = -sign(pk_old(gamma_xh));
0155 pk_old(gamma_xh) = sign(pk_old(gamma_xh)).*W_old(gamma_xh);
0156 pk = pk_old;
0157 dk = 0*pk;
0158 xk_1 = xh_old;
0159 
0160 % Initial step setup
0161 idelta = gamma_xh(1); flag = 1;    
0162 
0163 % initialize delx
0164 in_delx = [];
0165 delx_mode = opts.delx_mode;
0166 rhs = -dW.*z_x-Atdyu;
0167 
0168 if norm(xh_old)==0 && length(gamma_xh) == 1    
0169     update_mode = 'init0';
0170 else    
0171     update_mode = 'init1';
0172 end
0173 update_delx;
0174 
0175 %% loop parameters
0176 done = 0;
0177 iter = 0;
0178 itr_history = [];
0179 error_table = [];
0180 if err_record
0181     error_table = [epsilon err_fun(xk_1) length(gamma_xh)];
0182 end
0183 
0184 while iter < maxiter
0185     iter = iter+1;
0186     % warning('off','MATLAB:divideByZero')
0187     
0188     %% Homotopy
0189     x_k = xk_1;
0190     
0191     % Update direction
0192     delx_vec = zeros(N,1);
0193     delx_vec(gamma_xh) = delx;
0194     
0195     if ~isempty(idelta) && (sign(delx_vec(idelta)) == sign(pk_old(idelta)) && abs(x_k(idelta)) == 0)
0196         delta = 0; flag = 0;
0197     else
0198         pk = pk_old;
0199         % dk = AtA*delx_vec;
0200         % dk_temp = A*delx_vec;
0201         dk_temp = mvprod(A,delx_vec,gamma_xh,0);
0202         dk = A'*dk_temp+Atdyu;
0203         
0204         %%%--- compute step size
0205         in = [];
0206         
0207         % Setting shrinkage_flag to zero shrinks new active constraint towards the
0208         % final value instantly if doing so doesn't disturb the active set
0209         in.shrinkage_flag = 2; % determines how to select the stepsize/support
0210         in.nonneg = nonneg; % imposes non-negativity constraint on the solution
0211         in.pk = pk; in.dk = dk;
0212         in.ak = Wk; in.bk = dW;
0213         in.gamma = gamma_xh; in.gamma_c = gamma_xc;
0214         in.delx_vec = delx_vec; in.x = xk_1; 
0215         out = compute_delta(in);
0216         delta = out.delta; idelta = out.idelta;
0217         flag = out.flag;
0218     end
0219     e0 = epsilon;
0220     epsilon = e0 + delta;
0221     
0222     if epsilon > 1
0223         delta_end = 1-e0;
0224         xk_1 = x_k + delta_end*delx_vec;
0225         pk_old = pk + delta_end*dk;
0226         Wk = W_new;
0227         pk_old(gamma_xh) = sign(pk_old(gamma_xh)).*Wk(gamma_xh);
0228         break;
0229     end
0230     
0231     xk_1 = x_k + delta*delx_vec;
0232     gamma_old = gamma_xh;
0233     
0234     itr_history = [itr_history; idelta delta flag];
0235     
0236     % update support
0237     update_supp;
0238 
0239     pk_old = pk+delta*dk;
0240     Wk_old = Wk;
0241     Wk = (1-epsilon)*W_old + epsilon*W_new;
0242     pk_old([gamma_xh; idelta]) = sign(pk_old([gamma_xh; idelta])).*Wk([gamma_xh; idelta]);
0243     
0244     % update delx
0245     z_x = -sign(pk_old);
0246     rhs = -dW.*z_x-Atdyu;
0247     update_mode = 'update';
0248     update_delx;
0249     
0250     %     AtAgx = (A(:,gamma_xh)'*A(:,gamma_xh));
0251     %     delx2 = AtAgx\rhs(gamma_xh); % -AtAgx\(dW(gamma_xh).*z_x);
0252     %     figure(112); plot([delx delx2]);
0253     %     if norm(delx-delx2) > 1e-5
0254     %         stp = 1;
0255     %     end
0256     
0257     
0258     % Check convergence criterion (this can be useful)...
0259     if early_terminate
0260         if length(gamma_xh) < M/2
0261             xhat = zeros(N,1);
0262             % xhat(gamma_xh) = AtAgx\(A(:,gamma_xh)'*y);
0263             switch delx_mode
0264                 case 'mil'
0265                     xhat(gamma_xh) = iAtA*(A(:,gamma_xh)'*y);
0266                 case {'qr','chol'}
0267                     xhat(gamma_xh) = R\(R'\(A(:,gamma_xh)'*y));
0268             end
0269             if norm(y-A*xhat) < tau
0270                 xk_1 = xhat;
0271                 break;
0272             end
0273         end
0274     end
0275     
0276     
0277     if err_record        
0278         error_table = [error_table; epsilon err_fun(xk_1) length(gamma_xh)];
0279     end
0280     
0281     %% debug
0282     constr_violation = nnz((abs(pk_old(gamma_xc))-Wk(gamma_xc))>1e-10);
0283     sign_violation = nnz((sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh)))>1);
0284     if constr_violation
0285         chk = gamma_xc((abs(pk_old(gamma_xc))-Wk(gamma_xc))>1e-10);
0286         stp = 1;
0287         fprintf('problem... with constraint violation -- %s \n', mfilename);
0288         fprintf('Refactorize the matrix... recompute delx \n');
0289         % some times it comes here due to bad conditioning of AtAgx.
0290         update_mode = 'init0';
0291         update_delx;
0292     end
0293     if sign_violation>=1
0294         chk = gamma_xh(sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh))>1);
0295         stp = 1;
0296         fprintf('problem... sign mismatch -- %s \n',mfilename);
0297         fprintf('Refactorize the matrix... recompute delx \n');
0298         update_mode = 'init0';
0299         update_delx;
0300     end
0301     
0302     %% print and plot
0303     if mod(iter-1,verbose) == 0 && verbose
0304         fprintf(['iter = %d, delta = %3.4g, idelta = %d, flag = %d.\n'], iter, delta, idelta, flag);
0305     end
0306     if mod(iter-1,plots) == 0 && plots
0307         fig1 = figure(1234);
0308         subplot(2,1,1)
0309         hold off
0310         plot(pk,'.r', 'MarkerSize',14);
0311         hold on;
0312         plot(pk_old, 'LineWidth',1);
0313         
0314         if flag == 1
0315             plot(idelta, pk_old(idelta),'or','MarkerSize',18,'LineWidth',2);
0316             text(idelta, pk_old(idelta)*1.1, ['Incoming \gamma = ',num2str(idelta)],'FontSize',14);
0317         else
0318             plot(idelta, pk_old(idelta),'ok','MarkerSize',18,'LineWidth',2);
0319             text(idelta, pk_old(idelta)*1.1, ['Outgoing \gamma = ',num2str(idelta)],'FontSize',14);
0320         end
0321         set(gca,'FontSize',16, 'XLim',[1 N] );
0322         title(sprintf('BPDN shrinkage constraints: N = %d, M = %d', N, M));
0323         plot(1:N, Wk,'--k','MarkerSize',12);
0324         plot(1:N, -Wk, '--k','MarkerSize',12);
0325         plot(1:N, Wk_old,'--m','MarkerSize',12);
0326         plot(1:N, -Wk_old, '--m','MarkerSize',12);
0327         
0328         figure(fig1);
0329         subplot(2,1,2)
0330         hold off
0331         plot(x_k,'.r','MarkerSize',14); hold on;
0332         plot(xk_1,'LineWidth',1);
0333         if flag == 0
0334             plot(idelta, 0,'ok', 'MarkerSize',18,'LineWidth',2);
0335         else
0336             plot(idelta, 0,'or', 'MarkerSize',18,'LineWidth',2);
0337         end
0338         set(gca,'FontSize',16,'XLim',[1 N]);
0339         title(['Solution estimate at \epsilon = ',num2str(epsilon), ', iter. = ', num2str(iter)]);
0340         
0341         if iter == 1 && verbose
0342             disp('  ');
0343             disp('Every frame in the figure corresponds to a critical point on the homotopy path.')
0344             disp('Circle represents an incoming element, star represents an outgoing element.');
0345             disp(' ');
0346             disp('Put pause somewhere in the code to see this. ');
0347             disp('For now press some key to continue...');
0348             pause
0349         end
0350     end
0351     
0352 end
0353 
0354 %% debiasing step?
0355 if debias
0356     x_out = zeros(N,1);
0357     switch delx_mode
0358         case 'mil'
0359             x_out(gamma_xh) = iAtA*(A(:,gamma_xh)'*y);
0360         case {'qr','chol'}
0361             x_out(gamma_xh) = R\(R'\(A(:,gamma_xh)'*y));
0362     end
0363 else
0364     x_out = xk_1;
0365 end
0366 
0367 %
0368 if err_record
0369     error_table = [error_table; epsilon err_fun(x_out) length(gamma_xh)];
0370 end
0371 total_iter = iter;
0372 total_time = cputime-t0;
0373 
0374 %% Output the results
0375 out = opts;
0376 out.x_out = x_out;
0377 out.gamma = gamma_xh; % find(abs(xk_1)>0);
0378 out.iter = total_iter;
0379 out.time = total_time;
0380 out.error_table = error_table;
0381 out.pk = pk_old;
0382 switch delx_mode
0383     case 'mil'
0384         out.iAtA = iAtA;
0385     case 'qr'
0386         out.Q = Q;
0387         out.R = R;
0388     case 'chol'
0389         out.R = R;
0390     case 'qrM'
0391         out.Q0 = Q0;
0392         out.R0 = R0;
0393 end
0394