Quadratic Programming with Many Linear Constraints

This example shows how well the quadprog "active-set" algorithm performs in the presence of many linear constraints, as compared to the default "interior-point-convex" algorithm. Furthermore, the Lagrange multipliers from the "active-set" algorithm are exactly zero at inactive constraints, which can be helpful when you are looking for active constraints.

Problem Description

Create a pseudorandom quadratic problem with N variables and 10*N linear inequality constraints. Specify N = 150.

rng default % For reproducibility
N = 150;
rng default
A = randn([10*N,N]);
b = 10*ones(size(A,1),1);
f = sqrt(N)*rand(N,1);
H = 18*eye(N) + randn(N);
H = H + H';

Check that the resulting quadratic matrix is convex.

ee = min(eig(H))

ee = 
3.6976

All of the eigenvalues are positive, so the quadratic form x'*H*x is convex.

Include no linear equality constraints or bounds.

Aeq = [];
beq = [];
lb = [];
ub = [];

Solve Problem Using Two Algorithms

Set options to use the quadprog "active-set" algorithm. This algorithm requires an initial point. Set the initial point x0 to be a zero vector of length N.

opts = optimoptions("quadprog",Algorithm="active-set");
x0 = zeros(N,1);

Time the solution.

tic
[xa,fvala,eflaga,outputa,lambdaa] = quadprog(H,f,A,b,Aeq,beq,lb,ub,x0,opts);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

toc

Elapsed time is 0.069158 seconds.

Compare the solution time to the time of the default "interior-point-convex" algorithm.

tic
[xi,fvali,eflagi,outputi,lambdai] = quadprog(H,f,A,b);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

toc

Elapsed time is 3.924934 seconds.

The "active-set" algorithm is much faster on problems with many linear constraints.

Examine Lagrange Multipliers

The "active-set" algorithm reports only a few nonzero entries in the Lagrange multiplier structure associated with the linear constraint matrix.

nnz(lambdaa.ineqlin)

ans = 
14

In contrast, the "interior-point-convex" algorithm returns a Lagrange multiplier structure with all nonzero elements.

nnz(lambdai.ineqlin)

ans = 
1500

Nearly all of these Lagrange multipliers are smaller than N*eps in size.

nnz(abs(lambdai.ineqlin) > N*eps)

ans = 
23

In other words, the "active-set" algorithm gives clear indications of active constraints in the Lagrange multiplier structure, whereas the "interior-point-convex" algorithm does not.