nikcleju@2: % l1qc_logbarrier.m
nikcleju@2: %
nikcleju@2: % Solve quadratically constrained l1 minimization:
nikcleju@2: % min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
nikcleju@2: %
nikcleju@2: % Reformulate as the second-order cone program
nikcleju@2: % min_{x,u}  sum(u)   s.t.    x - u <= 0,
nikcleju@2: %                            -x - u <= 0,
nikcleju@2: %      1/2(||Ax-b||^2 - \epsilon^2) <= 0
nikcleju@2: % and use a log barrier algorithm.
nikcleju@2: %
nikcleju@2: % Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
nikcleju@2: %
nikcleju@2: % x0 - Nx1 vector, initial point.
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: % lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
nikcleju@2: %         Also, the number of log barrier iterations is completely
nikcleju@2: %         determined by lbtol.
nikcleju@2: %         Default = 1e-3.
nikcleju@2: %
nikcleju@2: % mu - Factor by which to increase the barrier constant at each iteration.
nikcleju@2: %      Default = 10.
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: function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)  
nikcleju@2: 
nikcleju@2: largescale = isa(A,'function_handle');
nikcleju@2: 
nikcleju@2: if (nargin < 6), lbtol = 1e-3; end
nikcleju@2: if (nargin < 7), mu = 10; end
nikcleju@2: if (nargin < 8), cgtol = 1e-8; end
nikcleju@2: if (nargin < 9), cgmaxiter = 200; end
nikcleju@2: 
nikcleju@2: newtontol = lbtol;
nikcleju@2: newtonmaxiter = 50;
nikcleju@2: 
nikcleju@2: N = length(x0);
nikcleju@2: 
nikcleju@2: % starting point --- make sure that it is feasible
nikcleju@2: if (largescale)
nikcleju@2:   if (norm(A(x0)-b) > epsilon)
nikcleju@2:     disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
nikcleju@2:     AAt = @(z) A(At(z));
nikcleju@2:     w = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
nikcleju@2:     if (cgres > 1/2)
nikcleju@2:       disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@2:       xp = x0;
nikcleju@2:       return;
nikcleju@2:     end
nikcleju@2:     x0 = At(w);
nikcleju@2:   end
nikcleju@2: else
nikcleju@2:   if (norm(A*x0-b) > epsilon)
nikcleju@2:     disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
nikcleju@2:     opts.POSDEF = true; opts.SYM = true;
nikcleju@2:     [w, hcond] = linsolve(A*A', b, opts);
nikcleju@2:     if (hcond < 1e-14)
nikcleju@2:       disp('A*At is ill-conditioned: cannot find starting point');
nikcleju@2:       xp = x0;
nikcleju@2:       return;
nikcleju@2:     end
nikcleju@2:     x0 = A'*w;
nikcleju@2:   end  
nikcleju@2: end
nikcleju@2: x = x0;
nikcleju@2: u = (0.95)*abs(x0) + (0.10)*max(abs(x0));
nikcleju@2: 
nikcleju@2: disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u)));
nikcleju@2: 
nikcleju@2: % choose initial value of tau so that the duality gap after the first
nikcleju@2: % step will be about the origial norm
nikcleju@2: tau = max((2*N+1)/sum(abs(x0)), 1);
nikcleju@2:                                                                                                                           
nikcleju@2: lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu));
nikcleju@2: disp(sprintf('Number of log barrier iterations = %d\n', lbiter));
nikcleju@2: 
nikcleju@2: totaliter = 0;
nikcleju@2: 
nikcleju@2: % Added by Nic
nikcleju@2: if lbiter == 0
nikcleju@2:     xp = zeros(size(x0));
nikcleju@2: end
nikcleju@2: 
nikcleju@2: for ii = 1:lbiter
nikcleju@2: 
nikcleju@2:   [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
nikcleju@2:   totaliter = totaliter + ntiter;
nikcleju@2:   
nikcleju@2:   disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
nikcleju@2:     ii, sum(abs(xp)), sum(up), tau, totaliter));
nikcleju@2:   
nikcleju@2:   x = xp;
nikcleju@2:   u = up;
nikcleju@2:  
nikcleju@2:   tau = mu*tau;
nikcleju@2:   
nikcleju@2: end
nikcleju@2: