% FORBES_QP
%
% FORBES_QP(H, q, A, lb, ub, Aeq, beq, lx, ux, opt, out1) solves the
% quadratic programming problem
%
% minimize (1/2)*x'*H*x + q'*x
% subject to Aeq*x == beq
% lb <= A*x <= ub
% lx <= x <= ux
%
% All arguments are optional and may be empty, except for the first.
% If last argument out1 is specified, then the solution process is
% warm-started based on the output of a previous call.
function out = forbes_qp(H, q, A, lb, ub, Aeq, beq, lx, ux, opt, out1)
t0 = tic();
% Arguments parsing
if nargin < 1 || isempty(H)
error('first parameter H is mandatory');
end
n = size(H, 2);
if nargin < 2 || isempty(q)
q = zeros(n, 1);
end
if nargin < 3 || isempty(A)
flag_ineq = 0;
else
flag_ineq = 1;
if size(A, 2) ~= n
error('argument A is incompatible with H and q');
end
end
m = size(A, 1);
if nargin < 4 || isempty(lb)
lb = -inf(m, 1);
end
if any(size(lb) ~= [m, 1])
error('argument lb is incompatible with A');
end
if nargin < 5 || isempty(ub)
ub = +inf(m, 1);
end
if any(size(ub) ~= [m, 1])
error('argument ub is incompatible with A');
end
if nargin < 6 || isempty(Aeq)
flag_eq = 0;
else
flag_eq = 1;
if size(Aeq, 2) ~= n
error('size of Aeq is incompatible with H and q');
end
end
if flag_eq == 1 && (nargin < 7 || isempty(beq) || any(size(beq) ~= [size(Aeq,1), 1]))
error('argument beq is incompatible with Aeq');
end
if nargin < 8 || isempty(lx)
flag_lx = 0;
lx = -inf(n, 1);
else
flag_lx = 1;
if isscalar(lx)
lx = lx*ones(n, 1);
end
end
if nargin < 9 || isempty(ux)
flag_ux = 0;
ux = +inf(n, 1);
else
flag_ux = 1;
if isscalar(ux)
ux = ux*ones(n, 1);
end
end
if nargin < 10, opt = []; end
if nargin < 11, out1 = []; end
if ~isfield(opt, 'prescale') || isempty(opt.prescale)
opt.prescale = true;
end
% Problem setup and solution
if flag_ineq == 0 && flag_eq == 0
f = quadratic(H, q);
g = indBox(lx, ux);
if isempty(out1)
x0 = zeros(n, 1);
else
opt.Lf = out1.solver.prob.Lf;
x0 = out1.x;
end
tprep = toc(t0);
out = forbes(f, g, x0, [], [], opt);
else
A_ext = A;
lb_ext = lb;
ub_ext = ub;
m_ext = m;
% Extend inequality constraints so as to include bounds on x
% (if necessary)
if flag_lx || flag_ux
A_ext = [A_ext; speye(n)];
lb_ext = [lb_ext; lx];
ub_ext = [ub_ext; ux];
m_ext = m_ext + n;
end
if flag_eq == 0 % NO equality constraints
f = quadratic(H, q);
opt_eigs.issym = 1;
opt_eigs.tol = 1e-3;
if (n <= 500 && min(eig(H)) <= 1e-16) || ...
(n > 500 && eigs(H, 1, 'SM', opt_eigs) <= 1e-16)
out.status = 2;
out.msg = 'not strongly convex';
return;
end
if opt.prescale
% Scale inequality constraints
scale = 1./sqrt(diag(A_ext*(H\A_ext')));
A_ext = diag(sparse(scale))*A_ext;
lb_ext = scale.*lb_ext;
ub_ext = scale.*ub_ext;
end
else % YES equality constraints
f = quadraticOverAffine(Aeq, beq, H, q);
if opt.prescale
scale = zeros(size(A_ext, 1), 1);
callfconj = f.makefconj();
[~, p] = callfconj(zeros(size(A_ext, 2),1));
for i = 1:size(A_ext, 1)
[~, dgradi] = callfconj(A_ext(i, :)');
w = A_ext(i, :)*(dgradi-p);
if w >= 1e-14
scale(i) = 1/sqrt(w);
else
scale(i) = 1;
end
end
% Scale inequality constraints
A_ext = diag(sparse(scale))*A_ext;
lb_ext = scale.*lb_ext;
ub_ext = scale.*ub_ext;
end
end
g = indBox(lb_ext, ub_ext);
constr = {A_ext, -speye(m_ext), zeros(m_ext, 1)};
if isempty(out1)
% cold start
y0 = zeros(m_ext, 1);
else
% warm start
opt.Lf = out1.solver.dual.prob.Lf;
y0 = [out1.y_ineq; out1.y_bnd];
end
tprep = toc(t0);
out_forbes = forbes(f, g, y0, [], constr, opt);
end
ttot = toc(t0);
out.status = out_forbes.flag;
out.msg = out_forbes.message;
out.x = out_forbes.x1;
out.y_ineq = out_forbes.y(1:m);
out.y_bnd = out_forbes.y(m+1:end);
out.pobj = (out.x'*(H*out.x))/2 + q'*out.x;
out.dobj = -out_forbes.dual.objective(end); % dual is solved as minimization
out.iterations = out_forbes.iterations;
out.preprocess = tprep;
out.time = ttot;
out.solver = out_forbes;
end