Quadratic Programming with Many Linear Constraints - MATLAB & Simulink (original) (raw)
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([10N,N]); b = 10ones(size(A,1),1); f = sqrt(N)rand(N,1); H = 18eye(N) + randn(N); H = H + H';
Check that the resulting quadratic matrix is convex.
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.
Elapsed time is 0.042058 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.
Elapsed time is 2.305694 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.
In contrast, the 'interior-point-convex' algorithm returns a Lagrange multiplier structure with all nonzero elements.
Nearly all of these Lagrange multipliers are smaller than N*eps in size.
nnz(abs(lambdai.ineqlin) > N*eps)
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.