Linearly Constrained Optimization
Mathwrist White Paper Series

Let $\psi (𝐱):{ℝ}^n\to ℝ$ be a twice continuously differentiable function with gradient $𝐠(𝐱)$ and Hessian $𝐇(𝐱)$. We want to solve the following optimization problem subject to a set of general linear constraints and simple bound constraints,

As we discussed in our Linear Programming (LP) and Quadratic Programming (QP) white papers, we rely on an active-set based framework to handle general linear constraints and simple bounds. The formulation (1) of linearly constrained (LC) optimization problems differs from a LP or QP problem only in the form of objective functions.

NPL LC optimization solvers naturally extend the active-set framework with specific customization of computing descent directions using line search methods. Same as LP and QP solvers, the LC solvers directly accept constraint specification in formulation (1). For brevity and without loss of generality, the discusion hereinafter assumes a canonical formulation,

Note that our implementation doesn’t convert formulation (1) to the canonical form (2). Please refer to our LP and QP white papers for the details of active-set method, i.e. how simple bounds and upper bounded linear constraints are handled. For the details of line search method, users can refer to our unconstrained optimization (UC) white paper.

We will use the same notation as our LP and QP white papers. In particular, $𝒲$ and $𝒩$ denote the index set of working constraints and non-working constraints respectively. $ℰ$ denotes the subset of equality constraints in LC problem (2). By construction, $ℰ\in 𝒲$. $𝐘$ and $𝐙$ represent the range space and null space matrix wrt active constraints ${𝐀}_{𝒲}$.

The solution ${𝐱}^{*}$ of a LC problem satisfies the following necessary and sufficient conditions:

1. 1.

   ${𝐱}^{*}$ is feasible. ${𝐀}_{𝒲}{𝐱}^{*}={𝐛}_{𝒲},{𝐀}_{𝒩}{𝐱}^{*}\ge {𝐛}_{𝒩}$;

2. 2.

   $𝐠({𝐱}^{*})={𝐀}_{𝒲}^T\lambda \text{ or equivalently }{𝐙}^T𝐠({𝐱}^{*})=0$;

3. 3.

   ${𝐙}^T𝐇({𝐱}^{*})𝐙$ is positive semi-definite (necessary) or positive definite (sufficient).

4. 4.

   Lagrange multiplier ${\lambda }_i\ge 0$ (necessary) or ${\lambda }_i>0$ (sufficient), $\forall i\in 𝒲-ℰ$;

Once an initial feasible point ${𝐱}_0$ is found, the active-set method automatically ensures all subsequent points ${𝐱}_k$ are feasible as well. The second condition is the stationary condition, effectively it requires the reduced gradient,

The third condition requires the reduced Hessian

to be at least positive semi-definite, also called curvature condition. The stationary and curvature conditions together determine if we can find a descent feasible direction wrt the current work set $𝒲$.

The last condition is the Lagrange multiplier condition. If , one can show that we can delete ${𝐀}_i$ from the current working set and obtain a descent feasible direction wrt the updated working set.

Given an initial guess $𝐱$, we use the “crash start” procedure in [2] to find an initial feasible point ${𝐱}_0$ and meanwhile identify an initial working set of constraints. This is the phase-I feasibility stage. Next, the algorithm enters iterations to repeatedly move along a descent feasible direction $𝐩$, update the working set and test optimality conditions.

At the $k$-th iteration, if ${𝐱}_k$ is not a stationary point wrt the current working set, or if $\widetilde{𝐇}({𝐱}_k)$ doesn’t have positive curvature, a descent direction $𝐩$ is computed as a null space direction $𝐩={\mathrm{𝐙𝐩}}_z$ for some ${𝐩}_z$. By moving along $𝐩$ with a step length $\alpha$, we first hit some non-working constraint $i\in 𝒩$. Now equality holds at ${𝐀}_i(𝐱+\alpha 𝐩)={𝐛}_i$. Constraint $i$ is then added to the working set $𝒲$. We update ${𝐱}_{k+1}={𝐱}_k+\alpha 𝐩$.

If ${𝐱}_k$ is a stationary point with positive curvature, optimality conditions are tested. A working constraint that violates the Lagrange multiplier condition is moved from $𝒲$ to $𝒩$. At this iteration, only the working set is updated. We then start the $k+1$ iteration with ${𝐱}_{k+1}={𝐱}_k$. Whenever the working set changes, i.e. constraints are added to or removed from $𝒲$, $𝐘$ and $𝐙$ are efficiently updated by a low-rank matrix update scheme, instead of being refactored from scratch.

For any search direction $𝐩$ and step length $\widehat{\alpha }$, by the Taylor expansion of $\psi (𝐱)$ around ${𝐱}_k$,

, where

We always move in a null space direction, $𝐩={\mathrm{𝐙𝐩}}_z$. The above expansion is written as

At a non-stationary point, we are able to compute a descent direction such that the linear term . However, when the reduced Hessian has positive curvature, the qudaratic term in equation (3) will dominate if $\widehat{\alpha }$ is too large. We use a line search method to compute an acceptable step length $\widehat{\alpha }$, which behaves as the upper bound of the actual step length $\alpha$.

Mathwrist NPL offers two LC solver classes,

1. 1.

   class GS_ActiveSetNewton

2. 2.

   class GS_ActiveSetQuasiNewton

Both solvers are hierarchically derived from the active-set framework base class. They all use the line search method to compute step length $\widehat{\alpha }$ and differ from each other on how to compute the step direction $𝐩$.

A special case of the linearly constrained optimization problem is when there are only simple bounds imposed to unknown variable $𝐱$,

, which is also called a box-constrained optimization problem. The overall logic to solve box-constrained problems is still following an active-set based strategy. But in this situation the working set $𝒲$ just holds a subset of variables that are fixed at their boundary. The non-working set $𝒩$ holds free variables.

Now searching in a null space direction $𝐩$ is simply to search in the subspace of free variables. We use two separate solver classes GS_SimpleBoundNewton and GS_SimpleBoundQuasiNewton to handle this special case more efficiently. Once a descent search direction $𝐩$ is computed, we still use a line search method to compute the upper bound of step length $\widehat{\alpha }$.

The special case solver class GS_SimpleBoundNewton computes the null space direction component ${𝐩}_z$ same as the general case in equation (6). Here the reduced gradient and reduced Hessian are just the free-variable parts of the gradient $𝐠({𝐱}_k)$ and Hessian $𝐇({𝐱}_k)$. Modified Cholesky is used to solve the linear equation.

The special case solver class GS_SimpleBoundQuasiNewton computes ${𝐩}_z$ as a Quasi-Newton direction (7). Again, ${𝐁}_k$ is a Hessian approximation by BFGS update. The reduced Hessian approximation ${\widetilde{𝐁}}_k$ is simply the free variable part of ${𝐁}_k$.

All four LC solvers

• •

  GS_ActiveSetNewton

• •

  GS_ActiveSetQuasiNewton

• •

  GS_SimpleBoundNewton

• •

  GS_SimpleBoundQuasiNewton

are all indirectly derived from the IterativeMethod base class. They all share the max number iteration control and convergence tolerance control defined in IterativeMethod base class. Let ${ϵ}_c$ be the convergence tolerance. At a stationary point, we test Lagrange multiplier condition by ${\lambda }_i\ge -{ϵ}_c$.

In addition, all four solvers use a stationary tolerance control ${ϵ}_s$ to determine whether the current iteration is at a stationary point. If the max-norm of reduced gradient satisfies , or if , we consider ${𝐱}_k$ is at a stationary point wrt the current working set of constraints.