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