activeset_bc

Fast active set algorithm for strictly convex quadratic program with bound constraints.

Fast active set algorithm for strictly convex quadratic program with bound constraints.

[x,lam,ws,oc,cc]=activeset_bc(H,g,LB,UB,tol,maxiter),
solves the strictly convex quadratic program of the form

min 0.5x'*H*x+g'*x, s.t. LB<=x<=UB,

where H=H'>0.

[x,lam,ws,oc,cc]=activeset_bc(H,g,LB,UB), work with default parameters.

This is modification of the algorithm 16.3 from Nocedal & Wright , p. 472.
Input:
H - positive definite and symmetric hessian matrix (nxn).
g - vetor (nx1).
tol - tolerance for equality constrained subproblem (16.39), p. 468 in N&W, (Default: 10^(-12))
maxiter - max. number of main iterations, (Default: 100*length(x))
Output:
x - optimal solution.
f=0.5x'*H*x+g'*x - minimum function value
lam(k) - Lagrange multiplier corresponding to LB(k)<=x(k)<=UB(k), k=1,..,n.
ws - is an optimal active set. If ws(k)=-1, then x(k)=LB(k). If ws(k)=+1,
then x(k)=UB(k), otherwise ws(k)=0.
oc - is an optimality measure i.e. ||Lx||=||H*x+g-A'*lam||, where
the Lagrangian

L=0.5x'*H*x+g'*x-lam'*(A*x-b),

and A=[I;-I], b=[LB;-UB].

cc - is complementarity measure i.e.

cc = max |lam(i)*min(x(i)-LB(i),UB(i)-x(i))|
i=1,..,n

The activeset_bc, is usually (but not always), faster than quadprog/interior-point
To compare with quadprog/interior-point try to run the following code:

clear;close all
n=1000; % Problem size.
ac=1.0; % Mean number of active constraints ~ ac (ac=0.625 --> 10%).
z=randn(n,1);H=eye(n)/ac+(z*z');g=randn(n,1); % Generation of the random problem.
LB=-ones(n,1);UB=ones(n,1); % Bounds.
t0=cputime;
[x,f,lam,ws,oc,cc]=activeset_bc(H,g,LB,UB); % Run activeset_bc.
dt=1000*(cputime-t0);nac=sum(abs(ws));
options=optimset('display','off');
t0=cputime;
x_qp=quadprog(H,g,[],[],[],[],LB,UB,[],options); % Run quadprog.
dt_qp=1000*(cputime-t0);dx=norm(x-x_qp);
disp(' ') % Compare the results
disp(['Problem size: ',num2str(n)])
disp(['Number of active constraints: ',num2str(nac)])
disp(['Norm of the difference between quadprog and activeset_bc: ',num2str(dx)])
disp(['Execution time quadprog: ',num2str(dt_qp), ' [ms]'])
disp(['Execution time activeset_bc: ',num2str(dt), ' [ms]'])
disp(['Optimality conditions: ',num2str(max(oc,cc))])
disp(' ')

(C) 2021, by dr Piotr Bania, [e-mail].
http://home.agh.edu.pl/~pba/