Quadratic Programming
Mathwrist White Paper Series

At each iteration, the active set method solves a sub QP problem subject to only equality constraints ${𝐀}_{𝒲}{𝐱}_k={𝐛}_{𝒲}$. If $\widetilde{𝐇}$ is positive definite, we can reach the local optimal ${𝐱}_k^{*}$ at a unit step $𝐩$. This search step $𝐩$ is the unique solution of the KKT system,

Note that the second equation in the KKT system now is ${𝐀}_{𝒲}𝐩=0$ to retain equality with respect to $𝒲$. In other words, $𝐩$ is a null space (optimality) step.

The KKT system (5) simultaneously solves the Lagrange multipliers $\lambda$ at ${𝐱}_k^{*}$ for all working constraints. If an inequality constraint $i$ in the working set, $i\in 𝒲-ℰ$, has negative Lagrange multiplier , we can remove constraint $i$ from $𝒲$ and start a new iteration. It can be shown that the search direction $𝐩$ computed from the KKT system (5) using the updated working set is a feasible direction to the deleted constraint $i$.

When we move along a null space direction $𝐩={\mathrm{𝐙𝐩}}_Z$ for some ${𝐪}_Z$, the objective function change is

The ${𝐙}^T𝐠$ term is oftern called reduced gradient, denoted as $\widetilde{𝐠}$ here. The change is,

After solving $𝐩$ from KKT system (5) and moving along this search direction, we could hit a non-working constraint before a unit step. Let $j\in 𝒩$ be the first blocking constraint among all non-working constraints. Let be the step length that we have moved along $𝐩$ to hit $j$. In this case we update ${𝐱}_{k+1}={𝐱}_k+\alpha 𝐩$. At this point, equality ${𝐀}_j{𝐱}_{k+1}={𝐛}_j$ exactly holds. We therefore remove $j$ from $𝒩$ and add it to $𝒲$. With this updated feasible point and new working set, we start the next iteration.

On the other hand, if we are able to make a unit step move $\alpha =1$, the local optimal ${𝐱}^{*}$ is reached. We test the Lagrange multipliers for all constraint $i\in 𝒲-ℰ$. Constraint $i$ with the most negative Lagrange is removed from the current working set $𝒲$ and added to $𝒩$. We update ${𝐱}_{k+1}={𝐱}^{*}=𝐱+𝐩$ and start the next iteration with the updated working set.

NPL QP active set solver directly works on the general formulation (1). At a local optimal point ${𝐱}^{*}$, the Lagrange multiplier requirement is ${\lambda }_i\ge 0$ for lower bounded constraints and ${\lambda }_i\le 0$ for upper bounded constraints. When we move along a descent feasible direction $𝐩$, for an inequality constraint $j\in 𝒩$, if , we check if the lower bound blocks the move. If ${𝐀}_j𝐩>0$, we check if its upper bound blocks. The working set $𝒲$ not only tracks the constraint index, but also which bound has equality being satisfied.

NPL handles the mixture of general linear constraints and simple bounds in the same active set framework. When a variable is equal to one of its boundary values, it is considered as a fixed variable. A variable stays strictly inside the constrained range is a free variable. Let ${𝐱}_b$ and ${𝐱}_f$ denote the fixed and free variables respectively. Let ${𝐁}_{𝒲}$ and ${𝐅}_{𝒲}$ be the working constraint coefficient sub matrix corresponding to ${𝐱}_b$ and ${𝐱}_f$. The set of working constraints effectively is

At each sub QP iteration, conceptually we have a KKT system in block structure

, where we have arranged the Hessian matrix $𝐇$ and gradient $𝐠$ into fixed variable and free variable blocks. The search direction ${𝐩}_f$ now is just relevant to free variables. At a local optimal ${𝐱}^{*}$, we now need also test the Lagrange multipliers ${\lambda }_b$ for the simple bound constraints ${𝐥}_b\le {𝐱}_b\le {𝐮}_b$.

From KKT system (8), clearly ${𝐩}_f$ is a null space direction to ${𝐅}_{𝒲}$. We maintain QR factorization of ${𝐅}_{𝒲}^T$,

Initially, a full QR decomposition is performed. Subsequently at each iteration, both simple bounds and general constraints can join or leave the working set. The block structure in (8) therefore changes accordingly. It is cost prohibitive to recompute a full QR factorization whenever ${𝐅}_{𝒲}$ changes. An efficient QR update scheme is implemented based on [2], section 5.2.4.

Because we use Cholesky factorization $\widetilde{𝐇}={\mathrm{𝐋𝐋}}^T$ to solve equation (3), it is also cost prohibitive to refactor $𝐋$ whenever $\widetilde{𝐇}$ is modified due to working set changes. We developed an efficient updating scheme to adopt the working set changes and economically obtain the new Cholesky factor from the previous Cholesky factor.

Based on [3], we implemented the inertia control strategy in combination with an extended KKT system. In summary, inertia control ensures that $\widetilde{𝐇}$ has at most one negative eigen value at any point of time. At a local optimal point, when a contraint $a_{*}$ is to be deleted from the working set, inertia control identifies if the deletion operation will cause $\widetilde{𝐇}$ being singular or indefinite. If so, $a_{*}$ is marked as deletion-pending but still kept in the extended KKT system. A zero curvature or negative curvature direction $𝐩$ is then computed accordingly.

When we move along $𝐩$ and hit a blocking constraint $a_i$, if adding $i$ to $𝒲$ lets the iteration reach to a vertex point, then clearly that is a local optimal point. We can remove $a_{*}$ from the extended KKT. If it is not a vertex, by testing the independency of $a_i$ wrt ${𝐀}_{𝒲}$ and $a_{*}$, inertia control can identify if such a contraint addition operation restores positive definiteness of $\widetilde{𝐇}$. If so, again $a_{*}$ can be physically removed from the extended KKT system. Otherwise, the algorithm keeps searching in zero or negative curvature directions and adding blocking constraints. Ultimately, either $\widetilde{𝐇}$ becomes positive definite or the current working set is a vertex point.

Most valuably in practical sense, [3] introduced a set of efficient schemes that can economically obtain the solutions of (9), (10) and (11) by updating the solutions of these linear systems at previous iterations. These update schemes are implemented in Mathwrist NPL.

When ${𝐱}_k$ is not a stationary point, a null space direction $𝐩$ that satisfies the following linear system is a zero curvature direction. Hence, the quadratic term in (6) vanishes.

Then we just need flip the sign if necessary so that the first term in equation (6) is negative.

When $\widetilde{𝐇}$ is indefinite at a non-stationary point, it is possible to find a negative curvature direction such that the quadratic term in (6) is negative. Note that we can again flip the sign of ${𝐩}_Z$ if necessary to make the first term in (6) negative as well. Such a direction $𝐩$ can be solved by the linear system

If $\widetilde{𝐇}$ is indefinite but ${𝐱}_k$ is a stationary point wrt the current working set, gradient $𝐠$ necessarily is a linear combination of working constraints ${𝐀}_{𝒲}$. Multiplying ${𝐙}^T$ to the right hand side of the first equation in (10) will be 0. So we cannot use (10) to find a negative curvature direction.

It can be shown that by following an inertia control strategy [3], the only possibility of having indefinite $\widetilde{𝐇}$ at a stationary point is by deleting a constraint ${𝐀}_i$ from the most recent local optimal point with ${\lambda }_i=0$. This observation allows us to find a negative curvature direction by solving,

Our active set framework always uses the “crash start” procedure [2] to find an initial feasible point and an initial working set $𝒲$. Then a full QR factorication is performed to ${𝐅}_{𝒲}^T$ to obtain the range space $𝐘$ and null space $𝐙$. Please refer to our LP white paper for details of this procedure.

For QP problems, the inertia control strategy requires the initial reduced Hessian ${𝐙}^T\mathrm{𝐇𝐙}$ being positive definite. If that is the case, surely for convex QP, we carry out an initial Cholesky factorization $\widetilde{𝐇}={\mathrm{𝐋𝐋}}^T$. The active set method is then ready to start its first iteration.

Given a general Hessian $𝐇$, it is possible that this initial $\widetilde{𝐇}$ is indefinite. We then perform a partial Cholesky [5],

, where the permutation $𝐏$ effectively arranges the null space as such that ${𝐙}_{+}^T{\mathrm{𝐇𝐙}}_{+}={𝐋}_{11}{\mathrm{𝐃𝐋}}_{11}^T$ is positive definite. All columns of ${𝐙}_{-}$ are treated as artificial constraints and added to the initial working set. The active set method then starts the first iteration. At a local optimal, these artificial constraints are deleted before any real constraint.