nikcleju@2: % l1qc_newton.m nikcleju@2: % nikcleju@2: % Newton algorithm for log-barrier subproblems for l1 minimization nikcleju@2: % with quadratic constraints. nikcleju@2: % nikcleju@2: % Usage: nikcleju@2: % [xp,up,niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, nikcleju@2: % newtontol, newtonmaxiter, cgtol, cgmaxiter) nikcleju@2: % nikcleju@2: % x0,u0 - starting points nikcleju@2: % nikcleju@2: % A - Either a handle to a function that takes a N vector and returns a K nikcleju@2: % vector , or a KxN matrix. If A is a function handle, the algorithm nikcleju@2: % operates in "largescale" mode, solving the Newton systems via the nikcleju@2: % Conjugate Gradients algorithm. nikcleju@2: % nikcleju@2: % At - Handle to a function that takes a K vector and returns an N vector. nikcleju@2: % If A is a KxN matrix, At is ignored. nikcleju@2: % nikcleju@2: % b - Kx1 vector of observations. nikcleju@2: % nikcleju@2: % epsilon - scalar, constraint relaxation parameter nikcleju@2: % nikcleju@2: % tau - Log barrier parameter. nikcleju@2: % nikcleju@2: % newtontol - Terminate when the Newton decrement is <= newtontol. nikcleju@2: % Default = 1e-3. nikcleju@2: % nikcleju@2: % newtonmaxiter - Maximum number of iterations. nikcleju@2: % Default = 50. nikcleju@2: % nikcleju@2: % cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix. nikcleju@2: % Default = 1e-8. nikcleju@2: % nikcleju@2: % cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored nikcleju@2: % if A is a matrix. nikcleju@2: % Default = 200. nikcleju@2: % nikcleju@2: % Written by: Justin Romberg, Caltech nikcleju@2: % Email: jrom@acm.caltech.edu nikcleju@2: % Created: October 2005 nikcleju@2: % nikcleju@2: nikcleju@2: nikcleju@2: function [xp, up, niter] = l1qc_newton(x0, u0, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter) nikcleju@2: nikcleju@2: % check if the matrix A is implicit or explicit nikcleju@2: largescale = isa(A,'function_handle'); nikcleju@2: nikcleju@2: % line search parameters nikcleju@2: alpha = 0.01; nikcleju@2: beta = 0.5; nikcleju@2: nikcleju@2: if (~largescale), AtA = A'*A; end nikcleju@2: nikcleju@2: % initial point nikcleju@2: x = x0; nikcleju@2: u = u0; nikcleju@2: if (largescale), r = A(x) - b; else r = A*x - b; end nikcleju@2: fu1 = x - u; nikcleju@2: fu2 = -x - u; nikcleju@2: fe = 1/2*(r'*r - epsilon^2); nikcleju@2: f = sum(u) - (1/tau)*(sum(log(-fu1)) + sum(log(-fu2)) + log(-fe)); nikcleju@2: nikcleju@2: niter = 0; nikcleju@2: done = 0; nikcleju@2: while (~done) nikcleju@2: nikcleju@2: if (largescale), atr = At(r); else atr = A'*r; end nikcleju@2: nikcleju@2: ntgz = 1./fu1 - 1./fu2 + 1/fe*atr; nikcleju@2: ntgu = -tau - 1./fu1 - 1./fu2; nikcleju@2: gradf = -(1/tau)*[ntgz; ntgu]; nikcleju@2: nikcleju@2: sig11 = 1./fu1.^2 + 1./fu2.^2; nikcleju@2: sig12 = -1./fu1.^2 + 1./fu2.^2; nikcleju@2: sigx = sig11 - sig12.^2./sig11; nikcleju@2: nikcleju@2: w1p = ntgz - sig12./sig11.*ntgu; nikcleju@2: if (largescale) nikcleju@2: h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr; nikcleju@2: [dx, cgres, cgiter] = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0); nikcleju@2: if (cgres > 1/2) nikcleju@2: disp('Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)'); nikcleju@2: xp = x; up = u; nikcleju@2: return nikcleju@2: end nikcleju@2: Adx = A(dx); nikcleju@2: else nikcleju@2: H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; nikcleju@2: opts.POSDEF = true; opts.SYM = true; nikcleju@2: [dx,hcond] = linsolve(H11p, w1p, opts); nikcleju@2: if (hcond < 1e-14) nikcleju@2: disp('Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)'); nikcleju@2: xp = x; up = u; nikcleju@2: return nikcleju@2: end nikcleju@2: Adx = A*dx; nikcleju@2: end nikcleju@2: du = (1./sig11).*ntgu - (sig12./sig11).*dx; nikcleju@2: nikcleju@2: % minimum step size that stays in the interior nikcleju@2: ifu1 = find((dx-du) > 0); ifu2 = find((-dx-du) > 0); nikcleju@2: aqe = Adx'*Adx; bqe = 2*r'*Adx; cqe = r'*r - epsilon^2; nikcleju@2: smax = min(1,min([... nikcleju@2: -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... nikcleju@2: (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe) nikcleju@2: ])); nikcleju@2: s = (0.99)*smax; nikcleju@2: nikcleju@2: % backtracking line search nikcleju@2: suffdec = 0; nikcleju@2: backiter = 0; nikcleju@2: while (~suffdec) nikcleju@2: xp = x + s*dx; up = u + s*du; rp = r + s*Adx; nikcleju@2: fu1p = xp - up; fu2p = -xp - up; fep = 1/2*(rp'*rp - epsilon^2); nikcleju@2: fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep)); nikcleju@2: flin = f + alpha*s*(gradf'*[dx; du]); nikcleju@2: suffdec = (fp <= flin); nikcleju@2: s = beta*s; nikcleju@2: backiter = backiter + 1; nikcleju@2: if (backiter > 32) nikcleju@2: disp('Stuck on backtracking line search, returning previous iterate. (See Section 4 of notes for more information.)'); nikcleju@2: xp = x; up = u; nikcleju@2: return nikcleju@2: end nikcleju@2: end nikcleju@2: nikcleju@2: % set up for next iteration nikcleju@2: x = xp; u = up; r = rp; nikcleju@2: fu1 = fu1p; fu2 = fu2p; fe = fep; f = fp; nikcleju@2: nikcleju@2: lambda2 = -(gradf'*[dx; du]); nikcleju@2: stepsize = s*norm([dx; du]); nikcleju@2: niter = niter + 1; nikcleju@2: done = (lambda2/2 < newtontol) | (niter >= newtonmaxiter); nikcleju@2: nikcleju@2: disp(sprintf('Newton iter = %d, Functional = %8.3f, Newton decrement = %8.3f, Stepsize = %8.3e', ... nikcleju@2: niter, f, lambda2/2, stepsize)); nikcleju@2: if (largescale) nikcleju@2: disp(sprintf(' CG Res = %8.3e, CG Iter = %d', cgres, cgiter)); nikcleju@2: else nikcleju@2: disp(sprintf(' H11p condition number = %8.3e', hcond)); nikcleju@2: end nikcleju@2: nikcleju@2: end nikcleju@2: nikcleju@2: nikcleju@2: nikcleju@2: