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