function [val,primals,duals] = iqp(H,c,A,b,verbose) % [val,primals,duals] = iqp(H,c,A,b,verbose) % % A simple interior point algorithm for solving the quadratic program % maximize - x'Hx/2 + c'x st Ax+b >= 0 % when H is positive semidefinite. Use H=[] for a linear program. % Use verbose=[] to turn off messages. % Detection of infeasibility/unboundedness is not very good: we just % return NaN if it appears that the optimal solution to either the % primal or the dual is very large. % General robustness could also be improved. % Copyright 1997 Geoff Gordon % a few parameters epsilon = 1e-8; % accuracy requirement for solution toobig = 1e+8; % give up if solution gets too large backoff = .99995; % don't go quite all the way to boundary maxiter = 75; % stop after this many iterations gamma = .5; % how fast to try to decrease barrier decay = .5; % -but effective gamma can be 1-(1-gamma)/decay mingamma = 1e-8; % -but bounded below by mingamma canreorder = 1; % set to zero to prevent reordering [m,n] = size(A); % default args if (nargin < 5) verbose = 1; end if (H == []) H = sparse(n,n); end % symmetric part of H is all that matters H = (H + H') / 2; % reorder for sparsity if (canreorder & issparse(A) & issparse(H)) order = symmmd(H + A' * A); H = H(order,order); A = A(:,order); c = c(order); else order = []; end % simple stupid initialization primals = zeros(n,1); duals = ones(m,1); slacks = ones(m,1); complementarity = (slacks' * duals) / m; % do until complementarity is small iter = 0; pstep = 1; dstep = 1; step = 0; barrier = 100; while ((complementarity > epsilon) & (iter < maxiter)) % target barrier is a fraction of complementarity determined % by how recently we've run against the edge of feasibility % (exact code here is arbitrary, anything reasonable should work) step = decay * min(pstep, dstep) * (step + 1/decay); effgamma = max(mingamma, (1 - step * (1 - gamma))); barrier = min(barrier, complementarity * effgamma); % compute rhs for Newton equations invslacks = 1 ./ slacks; d = spdiags(duals .* invslacks, 0, m, m); rhs = A' * (barrier .* invslacks + duals - d * b) + c; % solve Newton equations -- matlab is smart enough to use % Cholesky decomposition here newprimals = (H + (A' * d * A)) \ rhs; % from newprimals, compute newslacks newslacks = A * newprimals + b; % from newprimals and newslacks, compute newduals newduals = barrier .* invslacks - d * (newslacks - slacks); % compute primal and dual stepsizes: longest possible w/o hitting zero pstep = 1 ./ max(max(slacks - newslacks, slacks) ./ slacks); dstep = 1 ./ max(max(duals - newduals, duals) ./ duals); if (min(newslacks) <= 0) pstep = backoff * pstep; end if (min(newduals) <= 0) dstep = backoff * dstep; end % do the update primals = primals + pstep * (newprimals - primals); slacks = slacks + pstep * (newslacks - slacks); duals = duals + dstep * (newduals - duals); complementarity = (slacks' * duals) / m; iter = iter + 1; % give up if numbers get too big if max(max(slacks),max(duals)) > toobig val = NaN; if (order ~= []) primals(order) = primals; end if (verbose ~= []) disp('Warning: program appears to be ill-posed'); end return; end % report progress if (verbose ~= []) fprintf('iter %d, pstep %g, dstep %g, target %g, actual %g\n', ... iter, pstep, dstep, barrier, complementarity); end end val = primals' * c - primals' * H * primals / 2; if (order ~= []) primals(order) = primals; end