spgl1.m

Spectral Projected L1 solver

          if printTau || subspace
             printf(' %s',[tauFlag subFlag]);
          end
       end
       printf('\n');
    end
    printTau = false;
    subspace = false;
    
    % Update history info
    xNorm1(iter+1) = options.primal_norm(x,weights);
    rNorm2(iter+1) = rNorm;
    lambda(iter+1) = gNorm;
    
    if stat, break; end % Act on exit conditions.
        
    %==================================================================
    % Iterations begin here.
    %==================================================================
    iter = iter + 1;
    xOld = x;  fOld = f;  gOld = g;  rOld = r;

    try
       %---------------------------------------------------------------
       % Projected gradient step and linesearch.
       %---------------------------------------------------------------
       [f,x,r,nLine,stepG,lnErr] = ...
           spgLineCurvy(x,gStep*g,max(lastFv),@Aprod,b,@project,tau);
       nLineTot = nLineTot + nLine;
       if lnErr
          %  Projected backtrack failed. Retry with feasible dir'n linesearch.
          x    = xOld;
          dx   = project(x - gStep*g, tau) - x;
          gtd  = g'*dx;
          [f,x,r,nLine,lnErr] = spgLine(f,x,dx,gtd,max(lastFv),@Aprod,b);
          nLineTot = nLineTot + nLine;
       end
       if lnErr
       %  Failed again.  Revert to previous iterates and damp max BB step.
          if maxLineErrors <= 0
             stat = EXIT_LINE_ERROR;
          else
             stepMax = stepMax / 10;
             printf(['W: Linesearch failed with error %i. '...
                     'Damping max BB scaling to %6.1e.\n'],lnErr,stepMax);
             maxLineErrors = maxLineErrors - 1;
          end
       end

       %---------------------------------------------------------------
       % Subspace minimization (only if active-set change is small).
       %---------------------------------------------------------------
       doSubspaceMin = false;
       if subspaceMin
          g = - Aprod(r,2);
          [nnzX,nnzG,nnzIdx,nnzDiff] = activeVars(x,g,nnzOld,options);
          if ~nnzDiff
              if nnzX == nnzG, itnMaxLSQR = 20;
              else             itnMaxLSQR = 5;
              end
              nnzIdx = abs(x) >= optTol; 
              doSubspaceMin = true;
          end
       end

       if doSubspaceMin
     
          % LSQR parameters
          damp       = 1e-5;
          aTol       = 1e-1;
          bTol       = 1e-1;
          conLim     = 1e12;
          showLSQR   = 0;
       
          ebar   = sign(x(nnzIdx));
          nebar  = length(ebar);
          Sprod  = @(y,mode)LSQRprod(@Aprod,nnzIdx,ebar,n,y,mode);
       
          [dxbar, istop, itnLSQR] = ...
             lsqr(m,nebar,Sprod,r,damp,aTol,bTol,conLim,itnMaxLSQR,showLSQR);
              
          itnTotLSQR = itnTotLSQR + itnLSQR;
       
          if istop ~= 4  % LSQR iterations successful. Take the subspace step.
             % Push dx back into full space:  dx = Z dx.
             dx = zeros(n,1);
             dx(nnzIdx) = dxbar - (1/nebar)*(ebar'*dxbar)*dxbar;

             % Find largest step to a change in sign.
             block1 = nnzIdx  &  x < 0  &  dx > +pivTol;
             block2 = nnzIdx  &  x > 0  &  dx < -pivTol;
             alpha1 = Inf; alpha2 = Inf;
             if any(block1), alpha1 = min(-x(block1) ./ dx(block1)); end
             if any(block2), alpha2 = min(-x(block2) ./ dx(block2)); end
             alpha = min([1  alpha1  alpha2]);
             ensure(alpha >= 0);
             ensure(ebar'*dx(nnzIdx) <= optTol);          
          
             % Update variables.
             x    = x + alpha*dx;
             r    = b - Aprod(x,1);
             f    = r'*r / 2;
             subspace = true;
          end
       end
       
       ensure(options.primal_norm(x,weights) <= tau+optTol);

       %---------------------------------------------------------------
       % Update gradient and compute new Barzilai-Borwein scaling.
       %---------------------------------------------------------------
       g    = - Aprod(r,2);
       s    = x - xOld;
       y    = g - gOld;
       sts  = s'*s;
       sty  = s'*y;
       if   sty <= 0,  gStep = stepMax;
       else            gStep = min( stepMax, max(stepMin, sts/sty) );
       end
       
    catch % Detect matrix-vector multiply limit error
       err = lasterror;
       if strcmp(err.identifier,'SPGL1:MaximumMatvec')
         stat = EXIT_MATVEC_LIMIT;
         iter = iter - 1;
         x = xOld;  f = fOld;  g = gOld;  r = rOld;
         break;
       else
         rethrow(err);
       end
    end

    %------------------------------------------------------------------
    % Update function history.
    %------------------------------------------------------------------
    if singleTau || f > sigma^2 / 2 % Don't update if superoptimal.
       lastFv(mod(iter,nPrevVals)+1) = f;
       if fBest > f
          fBest = f;
          xBest = x;
       end
    end

end % while 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Restore best solution (only if solving single problem).
if singleTau && f > fBest
   rNorm = sqrt(2*fBest);
   printf('\n Restoring best iterate to objective %13.7e\n',rNorm);
   x = xBest;
   r = b - Aprod(x,1);
   g =   - Aprod(r,2);
   gNorm = options.dual_norm(g,weights);
   rNorm = norm(r,  2);
end

% Final cleanup before exit.
info.tau         = tau;
info.rNorm       = rNorm;
info.rGap        = rGap;
info.gNorm       = gNorm;
info.rGap        = rGap;
info.stat        = stat;
info.iter        = iter;
info.nProdA      = nProdA;
info.nProdAt     = nProdAt;
info.nNewton     = nNewton;
info.timeProject = timeProject;
info.timeMatProd = timeMatProd;
info.itnLSQR     = itnTotLSQR;
info.options     = options;
info.timeTotal   = toc;

info.xNorm1      = xNorm1(1:iter);
info.rNorm2      = rNorm2(1:iter);
info.lambda      = lambda(1:iter);

% Print final output.
switch (stat)
   case EXIT_OPTIMAL
      printf('\n EXIT -- Optimal solution found\n')
   case EXIT_ITERATIONS
      printf('\n ERROR EXIT -- Too many iterations\n');
   case EXIT_ROOT_FOUND
      printf('\n EXIT -- Found a root\n');
   case {EXIT_BPSOL1_FOUND, EXIT_BPSOL2_FOUND}
      printf('\n EXIT -- Found a BP solution\n');
   case EXIT_LINE_ERROR
      printf('\n ERROR EXIT -- Linesearch error (%i)\n',lnErr);
   case EXIT_SUBOPTIMAL_BP
      printf('\n EXIT -- Found a suboptimal BP solution\n');
   case EXIT_MATVEC_LIMIT
      printf('\n EXIT -- Maximum matrix-vector operations reached\n');
   case EXIT_ACTIVE_SET
      printf('\n EXIT -- Found a possible active set\n');
   otherwise
      error('Unknown termination condition\n');
end
printf('\n');
printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
   'Products with A',nProdA,'','Total time   (secs)',info.timeTotal);
printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
   'Products with A''',nProdAt,'','Project time (secs)',timeProject);
printf(' %-20s:  %6i %6s %-20s:  %6.1f\n',...
   'Newton iterations',nNewton,'','Mat-vec time (secs)',timeMatProd);
printf(' %-20s:  %6i %6s %-20s:  %6i\n', ...
   'Line search its',nLineTot,'','Subspace iterations',itnTotLSQR);
printf('\n');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NESTED FUNCTIONS.  These share some vars with workspace above.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
function z = Aprod(x,mode)
   if (nProdA + nProdAt >= maxMatvec)
     error('SPGL1:MaximumMatvec','');
   end
     
   tStart = toc;
   if mode == 1
      nProdA = nProdA + 1;
      if   explicit, z = A*x;
      else           z = A(x,1);
      end
   elseif mode == 2
      nProdAt = nProdAt + 1;
      if   explicit, z = A'*x;
      else           z = A(x,2);
      end
   else
      error('Wrong mode!');
   end
   timeMatProd = timeMatProd + (toc - tStart);
end % function Aprod

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function printf(varargin)
  if logLevel > 0
     fprintf(fid,varargin{:});
  end
end % function printf

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function x = project(x, tau)
   tStart      = toc;

   x = options.project(x,weights,tau);
   
   timeProject = timeProject + (toc - tStart);
end % function project

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of nested functions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end % function spg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PRIVATE FUNCTIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [nnzX,nnzG,nnzIdx,nnzDiff] = activeVars(x,g,nnzIdx,options)
% Find the current active set.
% nnzX    is the number of nonzero x.
% nnzG    is the number of elements in nnzIdx.
% nnzIdx  is a vector of primal/dual indicators.
% nnzDiff is the no. of elements that changed in the support.
  xTol    = min(.1,10*options.optTol);
  gTol    = min(.1,10*options.optTol);
  gNorm   = options.dual_norm(g,options.weights);
  nnzOld  = nnzIdx;

  % Reduced costs for postive & negative parts of x.
  z1 = gNorm + g;
  z2 = gNorm - g;

  % Primal/dual based indicators.
  xPos    = x >  xTol  &  z1 < gTol; %g < gTol;%
  xNeg    = x < -xTol  &  z2 < gTol; %g > gTol;%
  nnzIdx  = xPos | xNeg;

  % Count is based on simple primal indicator.
  nnzX    = sum(abs(x) >= xTol);
  nnzG    = sum(nnzIdx);
  
  if isempty(nnzOld)
     nnzDiff = inf;
  else
     nnzDiff = sum(nnzIdx ~= nnzOld);
  end
  
end % function spgActiveVars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function z = LSQRprod(Aprod,nnzIdx,ebar,n,dx,mode)
% Matrix multiplication for subspace minimization.
% Only called by LSQR.
  nbar = length(ebar);
   if mode == 1
      y = zeros(n,1);
      y(nnzIdx) = dx - (1/nbar)*(ebar'*dx)*ebar; % y(nnzIdx) = Z*dx
      z = Aprod(y,1);                            % z = S Z dx
   else
      y = Aprod(dx,2);
      z = y(nnzIdx) - (1/nbar)*(ebar'*y(nnzIdx))*ebar;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fNew,xNew,rNew,iter,err] = spgLine(f,x,d,gtd,fMax,Aprod,b)
% Nonmonotone linesearch.

EXIT_CONVERGED  = 0;
EXIT_ITERATIONS = 1;
maxIts = 10;
step   = 1;
iter   = 0;
gamma  = 1e-4;
gtd    = -abs(gtd); % 03 Aug 07: If gtd is complex,
                    % then should be looking at -abs(gtd).
while 1

    % Evaluate trial point and function value.
    xNew = x + step*d;
    rNew = b - Aprod(xNew,1);
    fNew = rNew'*rNew / 2;

    % Check exit conditions.
    if fNew < fMax + gamma*step*gtd  % Sufficient descent condition.
       err = EXIT_CONVERGED;
       break
    elseif  iter >= maxIts           % Too many linesearch iterations.
       err = EXIT_ITERATIONS;
       break
    end
    
    % New linesearch iteration.
    iter = iter + 1;
    
    % Safeguarded quadratic interpolation.
    if step <= 0.1
       step  = step / 2;
    else
       tmp = (-gtd*step^2) / (2*(fNew-f-step*gtd));
       if tmp < 0.1 || tmp > 0.9*step || isnan(tmp)
          tmp = step / 2;
       end
       step = tmp;
    end
    
end % while 1

end % function spgLine

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fNew,xNew,rNew,iter,step,err] = ...
    spgLineCurvy(x,g,fMax,Aprod,b,project,tau)
% Projected backtracking linesearch.
% On entry,
% g  is the (possibly scaled) steepest descent direction.

EXIT_CONVERGED  = 0;
EXIT_ITERATIONS = 1;
EXIT_NODESCENT  = 2;
gamma  = 1e-4;
maxIts = 10;
step   =  1;
sNorm  =  0;
scale  =  1;      % Safeguard scaling.  (See below.)
nSafe  =  0;      % No. of safeguarding steps.
iter   =  0;
debug  =  false;  % Set to true to enable log.
n      =  length(x);

if debug
   fprintf(' %5s  %13s  %13s  %13s  %8s\n',...
           'LSits','fNew','step','gts','scale');  
end
   
while 1

    % Evaluate trial point and function value.
    xNew     = project(x - step*scale*g, tau);
    rNew     = b - Aprod(xNew,1);
    fNew     = rNew'*rNew / 2;
    s        = xNew - x;
    gts      = scale * g' * s;
    if gts >= 0 % Should we check real and complex parts individually?
       err = EXIT_NODESCENT;
       break
    end

    if debug
       fprintf(' LS %2i  %13.7e  %13.7e  %13.6e  %8.1e\n',...
               iter,fNew,step,gts,scale);
    end
    
    % 03 Aug 07: If gts is complex, then should be looking at -abs(gts).
    if fNew < fMax - gamma*step*abs(gts)  % Sufficient descent condition.
       err = EXIT_CONVERGED;
       break
    elseif iter >= maxIts                 % Too many linesearch iterations.
       err = EXIT_ITERATIONS;
       break
    end
    
    % New linesearch iteration.
    iter = iter + 1;
    step = step / 2;

    % Safeguard: If stepMax is huge, then even damped search
    % directions can give exactly the same point after projection.  If
    % we observe this in adjacent iterations, we drastically damp the
    % next search direction.
    % 31 May 07: Damp consecutive safeguarding steps.
    sNormOld  = sNorm;
    sNorm     = norm(s) / sqrt(n);
    %   if sNorm >= sNormOld
    if abs(sNorm - sNormOld) <= 1e-6 * sNorm
       gNorm = norm(g) / sqrt(n);
       scale = sNorm / gNorm / (2^nSafe);
       nSafe = nSafe + 1;
    end
    
end % while 1

end % function spgLineCurvy