l1homotopy

 l1homotopy.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  \|W x\|_1 + 1/2*\|Ax-y\|_2^2 + (1-epsilon)u'x,

   u is defined as u = -W*sign(xh_old)-A'*(A*xh_old-y) 
   and xh_old is an arbitrary warm-start vector 
   (zero vector if no warm-start is available). 

   The homotopy is solved by changing epsilon from 0 to 1. 
   

 Optimality conditions:

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

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

  opts - input structure

   xh_old  -- old signal estimate (if warm-start is not provided, then xh_old is set to zero)
   gamma   -- support of xh_old

   pk_old  -- A'(A*xh_old-y)+u
   u       -- -W*sign(xh_old)-A'(A*xh_old-y)
   W       -- weights for the L1 term... 

   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
 Web: http://users.ece.gatech.edu/~sasif/

 References: 
 1. Sparse recovery of streaming signals using L1-homotopy, 
       by M. Salman Asif and Justin Romberg 
 2. Dynamic compressive sensing: Sparse recovery of streaming signals and video, 
       by M. Salman Asif (Ph.D. thesis) 

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

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