spgl1

SPGL1  Solve basis pursuit, basis pursuit denoise, and LASSO

 [x, r, g, info] = spgl1(A, b, tau, sigma, x0, options)

 ---------------------------------------------------------------------
 Solve the basis pursuit denoise (BPDN) problem

 (BPDN)   minimize  ||x||_1  subj to  ||Ax-b||_2 <= sigma,

 or the l1-regularized least-squares problem

 (LASSO)  minimize  ||Ax-b||_2  subj to  ||x||_1 <= tau.
 ---------------------------------------------------------------------

 INPUTS
 ======
 A        is an m-by-n matrix, explicit or an operator.
          If A is a function, then it must have the signature

          y = A(x,mode)   if mode == 1 then y = A x  (y is m-by-1);
                          if mode == 2 then y = A'x  (y is n-by-1).

 b        is an m-vector.
 tau      is a nonnegative scalar; see (LASSO).
 sigma    if sigma != inf or != [], then spgl1 will launch into a
          root-finding mode to find the tau above that solves (BPDN).
          In this case, it's STRONGLY recommended that tau = 0.
 x0       is an n-vector estimate of the solution (possibly all
          zeros). If x0 = [], then SPGL1 determines the length n via
          n = length( A'b ) and sets  x0 = zeros(n,1).
 options  is a structure of options from spgSetParms. Any unset options
          are set to their default value; set options=[] to use all
          default values.

 OUTPUTS
 =======
 x        is a solution of the problem
 r        is the residual, r = b - Ax
 g        is the gradient, g = -A'r
 info     is a structure with the following information:
          .tau     final value of tau (see sigma above)
          .rNorm   two-norm of the optimal residual
          .rGap    relative duality gap (an optimality measure)
          .gNorm   Lagrange multiplier of (LASSO)
          .stat    = 1 found a BPDN solution
                   = 2 found a BP sol'n; exit based on small gradient
                   = 3 found a BP sol'n; exit based on small residual
                   = 4 found a LASSO solution
                   = 5 error: too many iterations
                   = 6 error: linesearch failed
                   = 7 error: found suboptimal BP solution
                   = 8 error: too many matrix-vector products
          .time    total solution time (seconds)
          .nProdA  number of multiplications with A
          .nProdAt number of multiplications with A'

 OPTIONS
 =======
 Use the options structure to control various aspects of the algorithm:

 options.fid         File ID to direct log output
        .verbosity   0=quiet, 1=some output, 2=more output.
        .iterations  Max. number of iterations (default if 10*m).
        .bpTol       Tolerance for identifying a basis pursuit solution.
        .optTol      Optimality tolerance (default is 1e-4).
        .decTol      Larger decTol means more frequent Newton updates.
        .subspaceMin 0=no subspace minimization, 1=subspace minimization.

 EXAMPLE
 =======
   m = 120; n = 512; k = 20; % m rows, n cols, k nonzeros.
   p = randperm(n); x0 = zeros(n,1); x0(p(1:k)) = sign(randn(k,1));
   A  = randn(m,n); [Q,R] = qr(A',0);  A = Q';
   b  = A*x0 + 0.005 * randn(m,1);
   opts = spgSetParms('optTol',1e-4);
   [x,r,g,info] = spgl1(A, b, 0, 1e-3, [], opts); % Find BP sol'n.

 AUTHORS
 =======
  Ewout van den Berg (ewout78@cs.ubc.ca)
  Michael P. Friedlander (mpf@cs.ubc.ca)
    Scientific Computing Laboratory (SCL)
    University of British Columbia, Canada.

 BUGS
 ====
 Please send bug reports or comments to
            Michael P. Friedlander (mpf@cs.ubc.ca)
            Ewout van den Berg (ewout78@cs.ubc.ca)

0001 function [x,r,g,info] = spgl1( A, b, tau, sigma, x, options )
0002 %SPGL1  Solve basis pursuit, basis pursuit denoise, and LASSO
0003 %
0004 % [x, r, g, info] = spgl1(A, b, tau, sigma, x0, options)
0005 %
0006 % ---------------------------------------------------------------------
0007 % Solve the basis pursuit denoise (BPDN) problem
0008 %
0009 % (BPDN)   minimize  ||x||_1  subj to  ||Ax-b||_2 <= sigma,
0010 %
0011 % or the l1-regularized least-squares problem
0012 %
0013 % (LASSO)  minimize  ||Ax-b||_2  subj to  ||x||_1 <= tau.
0014 % ---------------------------------------------------------------------
0015 %
0016 % INPUTS
0017 % ======
0018 % A        is an m-by-n matrix, explicit or an operator.
0019 %          If A is a function, then it must have the signature
0020 %
0021 %          y = A(x,mode)   if mode == 1 then y = A x  (y is m-by-1);
0022 %                          if mode == 2 then y = A'x  (y is n-by-1).
0023 %
0024 % b        is an m-vector.
0025 % tau      is a nonnegative scalar; see (LASSO).
0026 % sigma    if sigma != inf or != [], then spgl1 will launch into a
0027 %          root-finding mode to find the tau above that solves (BPDN).
0028 %          In this case, it's STRONGLY recommended that tau = 0.
0029 % x0       is an n-vector estimate of the solution (possibly all
0030 %          zeros). If x0 = [], then SPGL1 determines the length n via
0031 %          n = length( A'b ) and sets  x0 = zeros(n,1).
0032 % options  is a structure of options from spgSetParms. Any unset options
0033 %          are set to their default value; set options=[] to use all
0034 %          default values.
0035 %
0036 % OUTPUTS
0037 % =======
0038 % x        is a solution of the problem
0039 % r        is the residual, r = b - Ax
0040 % g        is the gradient, g = -A'r
0041 % info     is a structure with the following information:
0042 %          .tau     final value of tau (see sigma above)
0043 %          .rNorm   two-norm of the optimal residual
0044 %          .rGap    relative duality gap (an optimality measure)
0045 %          .gNorm   Lagrange multiplier of (LASSO)
0046 %          .stat    = 1 found a BPDN solution
0047 %                   = 2 found a BP sol'n; exit based on small gradient
0048 %                   = 3 found a BP sol'n; exit based on small residual
0049 %                   = 4 found a LASSO solution
0050 %                   = 5 error: too many iterations
0051 %                   = 6 error: linesearch failed
0052 %                   = 7 error: found suboptimal BP solution
0053 %                   = 8 error: too many matrix-vector products
0054 %          .time    total solution time (seconds)
0055 %          .nProdA  number of multiplications with A
0056 %          .nProdAt number of multiplications with A'
0057 %
0058 % OPTIONS
0059 % =======
0060 % Use the options structure to control various aspects of the algorithm:
0061 %
0062 % options.fid         File ID to direct log output
0063 %        .verbosity   0=quiet, 1=some output, 2=more output.
0064 %        .iterations  Max. number of iterations (default if 10*m).
0065 %        .bpTol       Tolerance for identifying a basis pursuit solution.
0066 %        .optTol      Optimality tolerance (default is 1e-4).
0067 %        .decTol      Larger decTol means more frequent Newton updates.
0068 %        .subspaceMin 0=no subspace minimization, 1=subspace minimization.
0069 %
0070 % EXAMPLE
0071 % =======
0072 %   m = 120; n = 512; k = 20; % m rows, n cols, k nonzeros.
0073 %   p = randperm(n); x0 = zeros(n,1); x0(p(1:k)) = sign(randn(k,1));
0074 %   A  = randn(m,n); [Q,R] = qr(A',0);  A = Q';
0075 %   b  = A*x0 + 0.005 * randn(m,1);
0076 %   opts = spgSetParms('optTol',1e-4);
0077 %   [x,r,g,info] = spgl1(A, b, 0, 1e-3, [], opts); % Find BP sol'n.
0078 %
0079 % AUTHORS
0080 % =======
0081 %  Ewout van den Berg (ewout78@cs.ubc.ca)
0082 %  Michael P. Friedlander (mpf@cs.ubc.ca)
0083 %    Scientific Computing Laboratory (SCL)
0084 %    University of British Columbia, Canada.
0085 %
0086 % BUGS
0087 % ====
0088 % Please send bug reports or comments to
0089 %            Michael P. Friedlander (mpf@cs.ubc.ca)
0090 %            Ewout van den Berg (ewout78@cs.ubc.ca)
0091 
0092 % 15 Apr 07: First version derived from spg.m.
0093 %            Michael P. Friedlander (mpf@cs.ubc.ca).
0094 %            Ewout van den Berg (ewout78@cs.ubc.ca).
0095 % 17 Apr 07: Added root-finding code.
0096 % 18 Apr 07: sigma was being compared to 1/2 r'r, rather than
0097 %            norm(r), as advertised.  Now immediately change sigma to
0098 %            (1/2)sigma^2, and changed log output accordingly.
0099 % 24 Apr 07: Added quadratic root-finding code as an option.
0100 % 24 Apr 07: Exit conditions need to guard against small ||r||
0101 %            (ie, a BP solution).  Added test1,test2,test3 below.
0102 % 15 May 07: Trigger to update tau is now based on relative difference
0103 %            in objective between consecutive iterations.
0104 % 15 Jul 07: Added code to allow a limited number of line-search
0105 %            errors.
0106 % 23 Feb 08: Fixed bug in one-norm projection using weights. Thanks
0107 %            to Xiangrui Meng for reporting this bug.
0108 % 26 May 08: The simple call spgl1(A,b) now solves (BPDN) with sigma=0.
0109     
0110 %   spgl1.m
0111 %   $Id: spgl1.m 1225 2009-01-30 20:36:31Z ewout78 $
0112 %
0113 %   ----------------------------------------------------------------------
0114 %   This file is part of SPGL1 (Spectral Projected-Gradient for L1).
0115 %
0116 %   Copyright (C) 2007 Ewout van den Berg and Michael P. Friedlander,
0117 %   Department of Computer Science, University of British Columbia, Canada.
0118 %   All rights reserved. E-mail: <{ewout78,mpf}@cs.ubc.ca>.
0119 %
0120 %   SPGL1 is free software; you can redistribute it and/or modify it
0121 %   under the terms of the GNU Lesser General Public License as
0122 %   published by the Free Software Foundation; either version 2.1 of the
0123 %   License, or (at your option) any later version.
0124 %
0125 %   SPGL1 is distributed in the hope that it will be useful, but WITHOUT
0126 %   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
0127 %   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
0128 %   Public License for more details.
0129 %
0130 %   You should have received a copy of the GNU Lesser General Public
0131 %   License along with SPGL1; if not, write to the Free Software
0132 %   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
0133 %   USA
0134 %   ----------------------------------------------------------------------
0135 REVISION = '$Revision: 1017 $';
0136 DATE     = '$Date: 2008-06-16 22:43:07 -0700 (Mon, 16 Jun 2008) $';
0137 REVISION = REVISION(11:end-1);
0138 DATE     = DATE(35:50);
0139 
0140 tic;              % Start your watches!
0141 m = length(b);
0142 
0143 %----------------------------------------------------------------------
0144 % Check arguments.
0145 %----------------------------------------------------------------------
0146 if ~exist('options','var'), options = []; end
0147 if ~exist('x','var'), x = []; end
0148 if ~exist('sigma','var'), sigma = []; end
0149 if ~exist('tau','var'), tau = []; end
0150 
0151 if nargin < 2 || isempty(b) || isempty(A)
0152    error('At least two arguments are required');
0153 elseif isempty(tau) && isempty(sigma)
0154    tau = 0;
0155    sigma = 0;
0156    singleTau = false;
0157 elseif isempty(sigma) % && ~isempty(tau)  <-- implied
0158    singleTau = true;
0159 else
0160    if isempty(tau)
0161       tau = 0;
0162    end
0163    singleTau = false;
0164 end
0165 
0166 %----------------------------------------------------------------------
0167 % Grab input options and set defaults where needed.
0168 %----------------------------------------------------------------------
0169 defaultopts = spgSetParms(...
0170 'fid'        ,      1 , ... % File ID for output
0171 'verbosity'  ,      2 , ... % Verbosity level
0172 'iterations' ,   10*m , ... % Max number of iterations
0173 'nPrevVals'  ,      3 , ... % Number previous func values for linesearch
0174 'bpTol'      ,  1e-06 , ... % Tolerance for basis pursuit solution 
0175 'optTol'     ,  1e-04 , ... % Optimality tolerance
0176 'decTol'     ,  1e-04 , ... % Req'd rel. change in primal obj. for Newton
0177 'stepMin'    ,  1e-16 , ... % Minimum spectral step
0178 'stepMax'    ,  1e+05 , ... % Maximum spectral step
0179 'rootMethod' ,      2 , ... % Root finding method: 2=quad,1=linear (not used).
0180 'activeSetIt',    Inf , ... % Exit with EXIT_ACTIVE_SET if nnz same for # its.
0181 'subspaceMin',      0 , ... % Use subspace minimization
0182 'iscomplex'  ,    NaN , ... % Flag set to indicate complex problem
0183 'maxMatvec'  ,    Inf , ... % Maximum matrix-vector multiplies allowed
0184 'weights'    ,      1 , ... % Weights W in ||Wx||_1
0185 'project'    , @NormL1_project , ...
0186 'primal_norm', @NormL1_primal  , ...
0187 'dual_norm'  , @NormL1_dual      ...
0188    );
0189 options = spgSetParms(defaultopts, options);
0190 
0191 fid           = options.fid;
0192 logLevel      = options.verbosity;
0193 maxIts        = options.iterations;
0194 nPrevVals     = options.nPrevVals;
0195 bpTol         = options.bpTol;
0196 optTol        = options.optTol;
0197 decTol        = options.decTol;
0198 stepMin       = options.stepMin;
0199 stepMax       = options.stepMax;
0200 activeSetIt   = options.activeSetIt;
0201 subspaceMin   = options.subspaceMin;
0202 maxMatvec     = max(3,options.maxMatvec);
0203 weights       = options.weights;
0204 
0205 maxLineErrors = 10;     % Maximum number of line-search failures.
0206 pivTol        = 1e-12;  % Threshold for significant Newton step.
0207 
0208 %----------------------------------------------------------------------
0209 % Initialize local variables.
0210 %----------------------------------------------------------------------
0211 iter          = 0;  itnTotLSQR = 0; % Total SPGL1 and LSQR iterations.
0212 nProdA        = 0;  nProdAt    = 0;
0213 lastFv        = -inf(nPrevVals,1);  % Last m function values.
0214 nLineTot      = 0;                  % Total no. of linesearch steps.
0215 printTau      = false;
0216 nNewton       = 0;
0217 bNorm         = norm(b,2);
0218 stat          = false;
0219 timeProject   = 0;
0220 timeMatProd   = 0;
0221 nnzIter       = 0;                  % No. of its with fixed pattern.
0222 nnzIdx        = [];                 % Active-set indicator.
0223 subspace      = false;              % Flag if did subspace min in current itn.
0224 stepG         = 1;                  % Step length for projected gradient.
0225 testUpdateTau = 0;                  % Previous step did not update tau
0226 
0227 % Determine initial x, vector length n, and see if problem is complex
0228 explicit = ~(isa(A,'function_handle'));
0229 if isempty(x)
0230    if isnumeric(A)
0231       n = size(A,2);
0232       realx = isreal(A) && isreal(b);
0233    else
0234       x = Aprod(b,2);
0235       n = length(x);
0236       realx = isreal(x) && isreal(b);
0237    end
0238    x = zeros(n,1);
0239 else
0240    n     = length(x);
0241    realx = isreal(x) && isreal(b);
0242 end
0243 if isnumeric(A), realx = realx && isreal(A); end;
0244 
0245 % Override options when options.iscomplex flag is set
0246 if (~isnan(options.iscomplex)), realx = (options.iscomplex == 0); end
0247 
0248 % Check if all weights (if any) are strictly positive. In previous
0249 % versions we also checked if the number of weights was equal to
0250 % n. In the case of multiple measurement vectors, this no longer
0251 % needs to apply, so the check was removed.
0252 if ~isempty(weights)
0253   if any(~isfinite(weights))
0254      error('Entries in options.weights must be finite');
0255   end
0256   if any(weights <= 0)
0257      error('Entries in options.weights must be strictly positive');
0258   end
0259 else
0260   weights = 1;
0261 end
0262 
0263 % Quick exit if sigma >= ||b||.  Set tau = 0 to short-circuit the loop.
0264 if bNorm <= sigma
0265    printf('W: sigma >= ||b||.  Exact solution is x = 0.\n');
0266    tau = 0;  singleTau = true;
0267 end 
0268   
0269 % Don't do subspace minimization if x is complex.
0270 if ~realx && subspaceMin
0271    printf('W: Subspace minimization disabled when variables are complex.\n');
0272    subspaceMin = false;
0273 end
0274 
0275 % Pre-allocate iteration info vectors
0276 xNorm1 = zeros(min(maxIts,10000),1);
0277 rNorm2 = zeros(min(maxIts,10000),1);
0278 lambda = zeros(min(maxIts,10000),1);
0279 
0280 % Exit conditions (constants).
0281 EXIT_ROOT_FOUND    = 1;
0282 EXIT_BPSOL1_FOUND  = 2;
0283 EXIT_BPSOL2_FOUND  = 3;
0284 EXIT_OPTIMAL       = 4;
0285 EXIT_ITERATIONS    = 5;
0286 EXIT_LINE_ERROR    = 6;
0287 EXIT_SUBOPTIMAL_BP = 7;
0288 EXIT_MATVEC_LIMIT  = 8;
0289 EXIT_ACTIVE_SET    = 9; % [sic]
0290 
0291 %----------------------------------------------------------------------
0292 % Log header.
0293 %----------------------------------------------------------------------
0294 printf('\n');
0295 printf(' %s\n',repmat('=',1,80));
0296 printf(' SPGL1  v.%s (%s)\n', REVISION, DATE);
0297 printf(' %s\n',repmat('=',1,80));
0298 printf(' %-22s: %8i %4s'   ,'No. rows'          ,m       ,'');
0299 printf(' %-22s: %8i\n'     ,'No. columns'       ,n          );
0300 printf(' %-22s: %8.2e %4s' ,'Initial tau'       ,tau     ,'');
0301 printf(' %-22s: %8.2e\n'   ,'Two-norm of b'     ,bNorm      );
0302 printf(' %-22s: %8.2e %4s' ,'Optimality tol'    ,optTol  ,'');
0303 if singleTau
0304    printf(' %-22s: %8.2e\n'  ,'Target one-norm of x'  ,tau       );
0305 else
0306    printf(' %-22s: %8.2e\n','Target objective'  ,sigma      );
0307 end
0308 printf(' %-22s: %8.2e %4s' ,'Basis pursuit tol' ,bpTol   ,'');
0309 printf(' %-22s: %8i\n'     ,'Maximum iterations',maxIts     );
0310 printf('\n');
0311 if singleTau
0312    logB = ' %5i  %13.7e  %13.7e  %9.2e  %6.1f  %6i  %6i';
0313    logH = ' %5s  %13s  %13s  %9s  %6s  %6s  %6s\n';
0314    printf(logH,'Iter','Objective','Relative Gap','gNorm','stepG','nnzX','nnzG');
0315 else
0316    logB = ' %5i  %13.7e  %13.7e  %9.2e  %9.3e  %6.1f  %6i  %6i';
0317    logH = ' %5s  %13s  %13s  %9s  %9s  %6s  %6s  %6s  %13s\n';
0318    printf(logH,'Iter','Objective','Relative Gap','Rel Error',...
0319           'gNorm','stepG','nnzX','nnzG','tau');
0320 end    
0321     
0322 % Project the starting point and evaluate function and gradient.
0323 x         = project(x,tau);
0324 r         = b - Aprod(x,1);  % r = b - Ax
0325 g         =   - Aprod(r,2);  % g = -A'r
0326 f         = r'*r / 2; 
0327 
0328 % Required for nonmonotone strategy.
0329 lastFv(1) = f;
0330 fBest     = f;
0331 xBest     = x;
0332 fOld      = f;
0333 
0334 % Compute projected gradient direction and initial steplength.
0335 dx     = project(x - g, tau) - x;
0336 dxNorm = norm(dx,inf);
0337 if dxNorm < (1 / stepMax)
0338    gStep = stepMax;
0339 else
0340    gStep = min( stepMax, max(stepMin, 1/dxNorm) );
0341 end
0342 
0343 %----------------------------------------------------------------------
0344 % MAIN LOOP.
0345 %----------------------------------------------------------------------
0346 while 1
0347     
0348     %------------------------------------------------------------------
0349     % Test exit conditions.
0350     %------------------------------------------------------------------
0351 
0352     % Compute quantities needed for log and exit conditions.
0353     gNorm   = options.dual_norm(-g,weights);
0354     rNorm   = norm(r, 2);
0355     gap     = r'*(r-b) + tau*gNorm;
0356     rGap    = abs(gap) / max(1,f);
0357     aError1 = rNorm - sigma;
0358     aError2 = f - sigma^2 / 2;
0359     rError1 = abs(aError1) / max(1,rNorm);
0360     rError2 = abs(aError2) / max(1,f);
0361 
0362     % Count number of consecutive iterations with identical support.
0363     nnzOld = nnzIdx;
0364     [nnzX,nnzG,nnzIdx,nnzDiff] = activeVars(x,g,nnzIdx,options);
0365 
0366     if nnzDiff
0367        nnzIter = 0;
0368     else
0369        nnzIter = nnzIter + 1;
0370        if nnzIter >= activeSetIt, stat=EXIT_ACTIVE_SET; end
0371     end
0372     
0373     % Single tau: Check if we're optimal.
0374     % The 2nd condition is there to guard against large tau.
0375     if singleTau
0376        if rGap <= optTol || rNorm < optTol*bNorm
0377           stat  = EXIT_OPTIMAL;
0378        end
0379  
0380     % Multiple tau: Check if found root and/or if tau needs updating.
0381     else
0382        
0383        if rGap <= max(optTol, rError2) || rError1 <= optTol
0384           % The problem is nearly optimal for the current tau.
0385           % Check optimality of the current root.
0386           test1 = rNorm       <=   bpTol * bNorm;
0387           test2 = gNorm       <=   bpTol * rNorm;
0388           test3 = rError1     <=  optTol;
0389           test4 = rNorm       <=  sigma;
0390           
0391           if test4, stat=EXIT_SUBOPTIMAL_BP;end % Found suboptimal BP sol.
0392           if test3, stat=EXIT_ROOT_FOUND;   end % Found approx root.
0393           if test2, stat=EXIT_BPSOL2_FOUND; end % Gradient zero -> BP sol.
0394           if test1, stat=EXIT_BPSOL1_FOUND; end % Resid minim'zd -> BP sol.
0395        end
0396 
0397        testRelChange1 = (abs(f - fOld) <= decTol * f);
0398        testRelChange2 = (abs(f - fOld) <= 1e-1 * f * (abs(rNorm - sigma)));
0399        testUpdateTau  = ((testRelChange1 && rNorm >  2 * sigma) || ...
0400                          (testRelChange2 && rNorm <= 2 * sigma)) && ...
0401                          ~stat && ~testUpdateTau;
0402        
0403        if testUpdateTau
0404           % Update tau.
0405           tauOld   = tau;
0406           tau      = max(0,tau + (rNorm * aError1) / gNorm);
0407           nNewton  = nNewton + 1;
0408           printTau = abs(tauOld - tau) >= 1e-6 * tau; % For log only.
0409           if tau < tauOld
0410              % The one-norm ball has decreased.  Need to make sure that the
0411              % next iterate if feasible, which we do by projecting it.
0412              x = project(x,tau);
0413           end
0414        end
0415     end
0416 
0417     % Too many its and not converged.
0418     if ~stat  &&  iter >= maxIts
0419         stat = EXIT_ITERATIONS;
0420     end
0421 
0422     %------------------------------------------------------------------
0423     % Print log, update history and act on exit conditions.
0424     %------------------------------------------------------------------
0425     if logLevel >= 2 || singleTau || printTau || iter == 0 || stat
0426        tauFlag = '              '; subFlag = '';
0427        if printTau, tauFlag = sprintf(' %13.7e',tau);   end
0428        if subspace, subFlag = sprintf(' S %2i',itnLSQR); end
0429        if singleTau
0430           printf(logB,iter,rNorm,rGap,gNorm,log10(stepG),nnzX,nnzG);
0431           if subspace
0432              printf('  %s',subFlag);
0433           end
0434        else
0435           printf(logB,iter,rNorm,rGap,rError1,gNorm,log10(stepG),nnzX,nnzG);
0436           if printTau || subspace
0437              printf(' %s',[tauFlag subFlag]);
0438           end
0439        end
0440        printf('\n');
0441     end
0442     printTau = false;
0443     subspace = false;
0444     
0445     % Update history info
0446     xNorm1(iter+1) = options.primal_norm(x,weights);
0447     rNorm2(iter+1) = rNorm;
0448     lambda(iter+1) = gNorm;
0449     
0450     if stat, break; end % Act on exit conditions.
0451         
0452     %==================================================================
0453     % Iterations begin here.
0454     %==================================================================
0455     iter = iter + 1;
0456     xOld = x;  fOld = f;  gOld = g;  rOld = r;
0457 
0458     try
0459        %---------------------------------------------------------------
0460        % Projected gradient step and linesearch.
0461        %---------------------------------------------------------------
0462        [f,x,r,nLine,stepG,lnErr] = ...
0463            spgLineCurvy(x,gStep*g,max(lastFv),@Aprod,b,@project,tau);
0464        nLineTot = nLineTot + nLine;
0465        if lnErr
0466           %  Projected backtrack failed. Retry with feasible dir'n linesearch.
0467           x    = xOld;
0468           dx   = project(x - gStep*g, tau) - x;
0469           gtd  = g'*dx;
0470           [f,x,r,nLine,lnErr] = spgLine(f,x,dx,gtd,max(lastFv),@Aprod,b);
0471           nLineTot = nLineTot + nLine;
0472        end
0473        if lnErr
0474        %  Failed again.  Revert to previous iterates and damp max BB step.
0475           if maxLineErrors <= 0
0476              stat = EXIT_LINE_ERROR;
0477           else
0478              stepMax = stepMax / 10;
0479              printf(['W: Linesearch failed with error %i. '...
0480                      'Damping max BB scaling to %6.1e.\n'],lnErr,stepMax);
0481              maxLineErrors = maxLineErrors - 1;
0482           end
0483        end
0484 
0485        %---------------------------------------------------------------
0486        % Subspace minimization (only if active-set change is small).
0487        %---------------------------------------------------------------
0488        doSubspaceMin = false;
0489        if subspaceMin
0490           g = - Aprod(r,2);
0491           [nnzX,nnzG,nnzIdx,nnzDiff] = activeVars(x,g,nnzOld,options);
0492           if ~nnzDiff
0493               if nnzX == nnzG, itnMaxLSQR = 20;
0494               else             itnMaxLSQR = 5;
0495               end
0496               nnzIdx = abs(x) >= optTol; 
0497               doSubspaceMin = true;
0498           end
0499        end
0500 
0501        if doSubspaceMin
0502      
0503           % LSQR parameters
0504           damp       = 1e-5;
0505           aTol       = 1e-1;
0506           bTol       = 1e-1;
0507           conLim     = 1e12;
0508           showLSQR   = 0;
0509        
0510           ebar   = sign(x(nnzIdx));
0511           nebar  = length(ebar);
0512           Sprod  = @(y,mode)LSQRprod(@Aprod,nnzIdx,ebar,n,y,mode);
0513        
0514           [dxbar, istop, itnLSQR] = ...
0515              lsqr(m,nebar,Sprod,r,damp,aTol,bTol,conLim,itnMaxLSQR,showLSQR);
0516               
0517           itnTotLSQR = itnTotLSQR + itnLSQR;
0518        
0519           if istop ~= 4  % LSQR iterations successful. Take the subspace step.
0520              % Push dx back into full space:  dx = Z dx.
0521              dx = zeros(n,1);
0522              dx(nnzIdx) = dxbar - (1/nebar)*(ebar'*dxbar)*dxbar;
0523 
0524              % Find largest step to a change in sign.
0525              block1 = nnzIdx  &  x < 0  &  dx > +pivTol;
0526              block2 = nnzIdx  &  x > 0  &  dx < -pivTol;
0527              alpha1 = Inf; alpha2 = Inf;
0528              if any(block1), alpha1 = min(-x(block1) ./ dx(block1)); end
0529              if any(block2), alpha2 = min(-x(block2) ./ dx(block2)); end
0530              alpha = min([1  alpha1  alpha2]);
0531              ensure(alpha >= 0);
0532              ensure(ebar'*dx(nnzIdx) <= optTol);          
0533           
0534              % Update variables.
0535              x    = x + alpha*dx;
0536              r    = b - Aprod(x,1);
0537              f    = r'*r / 2;
0538              subspace = true;
0539           end
0540        end
0541        
0542        ensure(options.primal_norm(x,weights) <= tau+optTol);
0543 
0544        %---------------------------------------------------------------
0545        % Update gradient and compute new Barzilai-Borwein scaling.
0546        %---------------------------------------------------------------
0547        g    = - Aprod(r,2);
0548        s    = x - xOld;
0549        y    = g - gOld;
0550        sts  = s'*s;
0551        sty  = s'*y;
0552        if   sty <= 0,  gStep = stepMax;
0553        else            gStep = min( stepMax, max(stepMin, sts/sty) );
0554        end
0555        
0556     catch % Detect matrix-vector multiply limit error
0557        err = lasterror;
0558        if strcmp(err.identifier,'SPGL1:MaximumMatvec')
0559          stat = EXIT_MATVEC_LIMIT;
0560          iter = iter - 1;
0561          x = xOld;  f = fOld;  g = gOld;  r = rOld;
0562          break;
0563        else
0564          rethrow(err);
0565        end
0566     end
0567 
0568     %------------------------------------------------------------------
0569     % Update function history.
0570     %------------------------------------------------------------------
0571     if singleTau || f > sigma^2 / 2 % Don't update if superoptimal.
0572        lastFv(mod(iter,nPrevVals)+1) = f;
0573        if fBest > f
0574           fBest = f;
0575           xBest = x;
0576        end
0577     end
0578 
0579 end % while 1
0580 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0581 
0582 % Restore best solution (only if solving single problem).
0583 if singleTau && f > fBest
0584    rNorm = sqrt(2*fBest);
0585    printf('\n Restoring best iterate to objective %13.7e\n',rNorm);
0586    x = xBest;
0587    r = b - Aprod(x,1);
0588    g =   - Aprod(r,2);
0589    gNorm = options.dual_norm(g,weights);
0590    rNorm = norm(r,  2);
0591 end
0592 
0593 % Final cleanup before exit.
0594 info.tau         = tau;
0595 info.rNorm       = rNorm;
0596 info.rGap        = rGap;
0597 info.gNorm       = gNorm;
0598 info.rGap        = rGap;
0599 info.stat        = stat;
0600 info.iter        = iter;
0601 info.nProdA      = nProdA;
0602 info.nProdAt     = nProdAt;
0603 info.nNewton     = nNewton;
0604 info.timeProject = timeProject;
0605 info.timeMatProd = timeMatProd;
0606 info.itnLSQR     = itnTotLSQR;
0607 info.options     = options;
0608 info.timeTotal   = toc;
0609 
0610 info.xNorm1      = xNorm1(1:iter);
0611 info.rNorm2      = rNorm2(1:iter);
0612 info.lambda      = lambda(1:iter);
0613 
0614 % Print final output.
0615 switch (stat)
0616    case EXIT_OPTIMAL
0617       printf('\n EXIT -- Optimal solution found\n')
0618    case EXIT_ITERATIONS
0619       printf('\n ERROR EXIT -- Too many iterations\n');
0620    case EXIT_ROOT_FOUND
0621       printf('\n EXIT -- Found a root\n');
0622    case {EXIT_BPSOL1_FOUND, EXIT_BPSOL2_FOUND}
0623       printf('\n EXIT -- Found a BP solution\n');
0624    case EXIT_LINE_ERROR
0625       printf('\n ERROR EXIT -- Linesearch error (%i)\n',lnErr);
0626    case EXIT_SUBOPTIMAL_BP
0627       printf('\n EXIT -- Found a suboptimal BP solution\n');
0628    case EXIT_MATVEC_LIMIT
0629       printf('\n EXIT -- Maximum matrix-vector operations reached\n');
0630    case EXIT_ACTIVE_SET
0631       printf('\n EXIT -- Found a possible active set\n');
0632    otherwise
0633       error('Unknown termination condition\n');
0634 end
0635 printf('\n');
0636 printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
0637    'Products with A',nProdA,'','Total time   (secs)',info.timeTotal);
0638 printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
0639    'Products with A''',nProdAt,'','Project time (secs)',timeProject);
0640 printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
0641    'Newton iterations',nNewton,'','Mat-vec time (secs)',timeMatProd);
0642 printf(' %-20s:  %6i %6s %-20s:  %6i\n', ...
0643    'Line search its',nLineTot,'','Subspace iterations',itnTotLSQR);
0644 printf('\n');
0645 
0646 
0647 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0648 % NESTED FUNCTIONS.  These share some vars with workspace above.
0649 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0650     
0651 function z = Aprod(x,mode)
0652    if (nProdA + nProdAt >= maxMatvec)
0653      error('SPGL1:MaximumMatvec','');
0654    end
0655      
0656    tStart = toc;
0657    if mode == 1
0658       nProdA = nProdA + 1;
0659       if   explicit, z = A*x;
0660       else           z = A(x,1);
0661       end
0662    elseif mode == 2
0663       nProdAt = nProdAt + 1;
0664       if   explicit, z = A'*x;
0665       else           z = A(x,2);
0666       end
0667    else
0668       error('Wrong mode!');
0669    end
0670    timeMatProd = timeMatProd + (toc - tStart);
0671 end % function Aprod
0672 
0673 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0674 
0675 function printf(varargin)
0676   if logLevel > 0
0677      fprintf(fid,varargin{:});
0678   end
0679 end % function printf
0680 
0681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0682 
0683 function x = project(x, tau)
0684    tStart      = toc;
0685 
0686    x = options.project(x,weights,tau);
0687    
0688    timeProject = timeProject + (toc - tStart);
0689 end % function project
0690 
0691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0692 % End of nested functions.
0693 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0694 
0695 end % function spg
0696 
0697 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0698 % PRIVATE FUNCTIONS.
0699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0700 
0701 function [nnzX,nnzG,nnzIdx,nnzDiff] = activeVars(x,g,nnzIdx,options)
0702 % Find the current active set.
0703 % nnzX    is the number of nonzero x.
0704 % nnzG    is the number of elements in nnzIdx.
0705 % nnzIdx  is a vector of primal/dual indicators.
0706 % nnzDiff is the no. of elements that changed in the support.
0707   xTol    = min(.1,10*options.optTol);
0708   gTol    = min(.1,10*options.optTol);
0709   gNorm   = options.dual_norm(g,options.weights);
0710   nnzOld  = nnzIdx;
0711 
0712   % Reduced costs for postive & negative parts of x.
0713   z1 = gNorm + g;
0714   z2 = gNorm - g;
0715 
0716   % Primal/dual based indicators.
0717   xPos    = x >  xTol  &  z1 < gTol; %g < gTol;%
0718   xNeg    = x < -xTol  &  z2 < gTol; %g > gTol;%
0719   nnzIdx  = xPos | xNeg;
0720 
0721   % Count is based on simple primal indicator.
0722   nnzX    = sum(abs(x) >= xTol);
0723   nnzG    = sum(nnzIdx);
0724   
0725   if isempty(nnzOld)
0726      nnzDiff = inf;
0727   else
0728      nnzDiff = sum(nnzIdx ~= nnzOld);
0729   end
0730   
0731 end % function spgActiveVars
0732 
0733 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0734 
0735 function z = LSQRprod(Aprod,nnzIdx,ebar,n,dx,mode)
0736 % Matrix multiplication for subspace minimization.
0737 % Only called by LSQR.
0738   nbar = length(ebar);
0739    if mode == 1
0740       y = zeros(n,1);
0741       y(nnzIdx) = dx - (1/nbar)*(ebar'*dx)*ebar; % y(nnzIdx) = Z*dx
0742       z = Aprod(y,1);                            % z = S Z dx
0743    else
0744       y = Aprod(dx,2);
0745       z = y(nnzIdx) - (1/nbar)*(ebar'*y(nnzIdx))*ebar;
0746    end
0747 end
0748 
0749 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0750 
0751 function [fNew,xNew,rNew,iter,err] = spgLine(f,x,d,gtd,fMax,Aprod,b)
0752 % Nonmonotone linesearch.
0753 
0754 EXIT_CONVERGED  = 0;
0755 EXIT_ITERATIONS = 1;
0756 maxIts = 10;
0757 step   = 1;
0758 iter   = 0;
0759 gamma  = 1e-4;
0760 gtd    = -abs(gtd); % 03 Aug 07: If gtd is complex,
0761                     % then should be looking at -abs(gtd).
0762 while 1
0763 
0764     % Evaluate trial point and function value.
0765     xNew = x + step*d;
0766     rNew = b - Aprod(xNew,1);
0767     fNew = rNew'*rNew / 2;
0768 
0769     % Check exit conditions.
0770     if fNew < fMax + gamma*step*gtd  % Sufficient descent condition.
0771        err = EXIT_CONVERGED;
0772        break
0773     elseif  iter >= maxIts           % Too many linesearch iterations.
0774        err = EXIT_ITERATIONS;
0775        break
0776     end
0777     
0778     % New linesearch iteration.
0779     iter = iter + 1;
0780     
0781     % Safeguarded quadratic interpolation.
0782     if step <= 0.1
0783        step  = step / 2;
0784     else
0785        tmp = (-gtd*step^2) / (2*(fNew-f-step*gtd));
0786        if tmp < 0.1 || tmp > 0.9*step || isnan(tmp)
0787           tmp = step / 2;
0788        end
0789        step = tmp;
0790     end
0791     
0792 end % while 1
0793 
0794 end % function spgLine
0795 
0796 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0797 
0798 function [fNew,xNew,rNew,iter,step,err] = ...
0799     spgLineCurvy(x,g,fMax,Aprod,b,project,tau)
0800 % Projected backtracking linesearch.
0801 % On entry,
0802 % g  is the (possibly scaled) steepest descent direction.
0803 
0804 EXIT_CONVERGED  = 0;
0805 EXIT_ITERATIONS = 1;
0806 EXIT_NODESCENT  = 2;
0807 gamma  = 1e-4;
0808 maxIts = 10;
0809 step   =  1;
0810 sNorm  =  0;
0811 scale  =  1;      % Safeguard scaling.  (See below.)
0812 nSafe  =  0;      % No. of safeguarding steps.
0813 iter   =  0;
0814 debug  =  false;  % Set to true to enable log.
0815 n      =  length(x);
0816 
0817 if debug
0818    fprintf(' %5s  %13s  %13s  %13s  %8s\n',...
0819            'LSits','fNew','step','gts','scale');  
0820 end
0821    
0822 while 1
0823 
0824     % Evaluate trial point and function value.
0825     xNew     = project(x - step*scale*g, tau);
0826     rNew     = b - Aprod(xNew,1);
0827     fNew     = rNew'*rNew / 2;
0828     s        = xNew - x;
0829     gts      = scale * g' * s;
0830     if gts >= 0 % Should we check real and complex parts individually?
0831        err = EXIT_NODESCENT;
0832        break
0833     end
0834 
0835     if debug
0836        fprintf(' LS %2i  %13.7e  %13.7e  %13.6e  %8.1e\n',...
0837                iter,fNew,step,gts,scale);
0838     end
0839     
0840     % 03 Aug 07: If gts is complex, then should be looking at -abs(gts).
0841     if fNew < fMax - gamma*step*abs(gts)  % Sufficient descent condition.
0842        err = EXIT_CONVERGED;
0843        break
0844     elseif iter >= maxIts                 % Too many linesearch iterations.
0845        err = EXIT_ITERATIONS;
0846        break
0847     end
0848     
0849     % New linesearch iteration.
0850     iter = iter + 1;
0851     step = step / 2;
0852 
0853     % Safeguard: If stepMax is huge, then even damped search
0854     % directions can give exactly the same point after projection.  If
0855     % we observe this in adjacent iterations, we drastically damp the
0856     % next search direction.
0857     % 31 May 07: Damp consecutive safeguarding steps.
0858     sNormOld  = sNorm;
0859     sNorm     = norm(s) / sqrt(n);
0860     %   if sNorm >= sNormOld
0861     if abs(sNorm - sNormOld) <= 1e-6 * sNorm
0862        gNorm = norm(g) / sqrt(n);
0863        scale = sNorm / gNorm / (2^nSafe);
0864        nSafe = nSafe + 1;
0865     end
0866     
0867 end % while 1
0868 
0869 end % function spgLineCurvy