l1qc_newton.m

基于内点法的解l2_l1和l2_TV优化问题的matlab代码

% l1qc_newton.m
%
% Newton algorithm for log-barrier subproblems for l1 minimization
% with quadratic constraints.
%
% Usage: 
% [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, 
%                             newtontol, newtonmaxiter, cgtol, cgmaxiter)
%
% x0,u0 - starting points
%
% A - Either a handle to a function that takes a N vector and returns a K 
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% epsilon - scalar, constraint relaxation parameter
%
% tau - Log barrier parameter.
%
% newtontol - Terminate when the Newton decrement is <= newtontol.
%         Default = 1e-3.
%
% newtonmaxiter - Maximum number of iterations.
%         Default = 50.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: [e-mail]
% Created: October 2005
%


function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) 

% check if the matrix A is implicit or explicit
largescale = isa(A,'function_handle');

% line search parameters
alpha = 0.01;
beta = 0.5;  

N = length(x0);

if (~largescale), AtA = A'*A; end

% initial point
x = x0;
u = u0;
if (largescale), r = A(x) - b; else, r = A*x - b; end
fu1 = x - u;
fu2 = -x - u;
fe = 1/2*(r'*r - epsilon^2);
f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe));

niter = 0;
done = 0;
while (~done)
  
  if (largescale), atr = At(r); else, atr = A'*r; end
  
  ntgz = 1./fu1 - 1./fu2 + 1/fe*atr;
  ntgu = -tau - 1./fu1 - 1./fu2;
  gradf = -(1/tau)*[ntgz; ntgu];
  
  sig11 = 1./fu1.^2 + 1./fu2.^2;
  sig12 = -1./fu1.^2 + 1./fu2.^2;
  sigx = sig11 - sig12.^2./sig11;
    
  w1p = ntgz - sig12./sig11.*ntgu;
  if (largescale)
    h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr;
    [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0);
    Adx = A(dx);
  else
    H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; 
    dx = H11p\w1p;
    Adx = A*dx;
  end
  du = (1./sig11).*ntgu - (sig12./sig11).*dx;  
 
  % minimum step size that stays in the interior
  s = 1;
  xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
  coneiter = 0;
  while ( (max(abs(xp)-up) > 0) | (rp'*rp > epsilon^2) )
    s = beta*s;
    xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
    coneiter = coneiter + 1;
    if (coneiter > 32)
      disp('Stuck on cone iterations, returning previous iterate.');
      xp = x;  up = u;
      return
    end    
  end
    
  % backtracking line search
  fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);
  fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));
  flin = f + alpha*s*(gradf'*[dx; du]);
  backiter = 0;
  while (fp > flin)
    s = beta*s;
    xp = x + s*dx;  up = u + s*du;  rp = r + s*Adx;
    fu1p = xp - up;  fu2p = -xp - up;  fep = 1/2*(rp'*rp - epsilon^2);
    fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep));
    flin = f + alpha*s*(gradf'*[dx; du]);
    backiter = backiter + 1;
    if (backiter > 32)
      disp('Stuck on backtracking line search, returning previous iterate.');
      xp = x;  up = u;
      return
    end
  end
  
  % set up for next iteration
  x = xp; u = up;  r = rp;
  fu1 = fu1p;  fu2 = fu2p;  fe = fep;  f = fp;
  
  lambda2 = -(gradf'*[dx; du]);
  stepsize = s*norm([dx; du]);
  niter = niter + 1;
  done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter);
  
  disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e, Cone iterations = %8.3f, Backtrack iterations = %d', ...
    niter, f, lambda2/2, stepsize, coneiter, backiter));
  if (largescale)
    disp(sprintf('                  CG Res = %8.3e, CG Iter = %d', cgres, cgiter));
  end
      
end