Linear Programming
Mathwrist White Paper Series

(January 20, 2023)

Abstract

This document presents a technical overview of the Linear Programming (LP) feature in Mathwrist’s C++ Numerical Programming Library (NPL). Two algorithms are implemented to solve LP problems. One is the classic Simplex method. The second one is an active-set based method with linear objective function.

1 Introduction

In canonical form, linear programming (LP) solves the following minimization problem with a linear objective function subject to linear constraints.

\min_{\mathbf{x}\in\mathbb{R}^{n}}\mathbf{c}^{T}\mathbf{x},\mbox{ s.t. }% \mathbf{A}\mathbf{x}=\mathbf{b}\mbox{ and }\mathbf{x}\geq 0

(1)

To better support practical use, we actually formulate LP problems in a more general form,

\min_{\mathbf{x}\in\mathbb{R}^{n}}\mathbf{c}^{T}\mathbf{x},\mbox{ s.t. }% \mathbf{b}_{l}\leq\mathbf{A}\mathbf{x}\leq\mathbf{b}_{u}\mbox{ and }\mathbf{l}% \leq\mathbf{x}\leq\mathbf{u}

(2)

We shall refer $\mathbf{c}$ as the cost vector, $\mathbf{b}_{l}\leq\mathbf{A}\mathbf{x}\leq\mathbf{b}_{u}$ as general linear constraints and $\mathbf{l}\leq\mathbf{x}\leq\mathbf{u}$ as simple bound (box) constraints respectively. Our implementation of both the Simplex method and the active-set method all take this general LP formulation at user interface level. Internally, the Simplex method needs convert this general form to the canonical form.

In the general form, let $\mathcal{E}$ denote the index subset of general constraints that require equality $\mathbf{A}_{\mathcal{E}}\mathbf{x}=\mathbf{b}_{\mathcal{E}}$ to hold. Assume $\mathbf{A}_{\mathcal{E}}$ has full row rank. Clearly, the number of equality constraints should be less than $n$ . Otherwise the solution is purely determined by $\mathbf{A}_{\mathcal{E}}$ .

Our main references for Simplex method are [1] by David G. Luenberger and [2] by Jorge Nocedal and Stephen J. Wright. Inspired by [3] by Gill et al., we effectively implemented an active-set based framework that can be customized to handle LP or many other problems.

2 Simplex Method

The classic Simplex algorithm is derived from the canonical LP formulation. In this canonical form, all general constraints are equality constraints. We shall assume $\mathbf{A}$ has full row rank. Let $m$ be the number of linear equality contraints and $n$ be the number of unknowns. Again, $m<n$ , otherwise we either don’t have a feasible solution, or the solution is just determined by the linear system $\mathbf{A}\mathbf{x}=\mathbf{b}$ .

Next, let $\mathcal{B}$ be an index set of exact $m$ number of index elements, $\mathcal{B}\in\{0,1,\cdots,n\}$ . Let $\mathcal{N}$ denote the complement subset. A unknown variable $x_{i}$ that has index $i\in\mathcal{B}$ or $i\in\mathcal{N}$ is called a basic variable or non-basic variable respectively. We can shuffle the order of unknown variables and rewrite the linear equality constraints as

\displaystyle\left(\begin{array}[]{c|c}\hidden@noalign{}\hfil\mathbf{B}&% \mathbf{N}\end{array}\right)\left(\begin{array}[]{c}\hidden@noalign{}\hfil% \mathbf{x}_{B}\cr\hidden@noalign{}\hfil\mathbf{x}_{N}\end{array}\right)=% \mathbf{b}

, where $\mathbf{x}_{B}$ and $\mathbf{x}_{N}$ are basic variables and non-basic variables. $\mathbf{B}$ and $\mathbf{N}$ are the constraint coefficient submatrix for $\mathbf{x}_{B}$ and $\mathbf{x}_{N}$ respectively. If we could find an index set $\mathcal{B}$ such that $\mathbf{B}$ is nonsingular and $\mathbf{x}_{B}=\mathbf{B}^{-1}\mathbf{b}\geq 0$ , we can then fix $\mathbf{x}_{N}=0$ to obtaint a basic feasible solution.

The Fundamental Theorem of Linear Programming, i.e. [2] section 13.2, ensures that if the LP problem (1) has solutions, at least one of them is a basic point. Given an initial basic feasible point $\mathbf{x}_{0}$ and the associated index set, Simplex method iteratively searches in the space of basic feasible points. At each iteration, it tests whether the optimality conditions are satisfied. If not, a non-basic index $q\in\mathcal{N}$ is identified as the “enter index” and a basic index $p\in\mathcal{B}$ is identified as the “quit index”. By increasing variable $x_{q}$ from zero and decreasing variable $x_{p}$ to 0, the objective function $\mathbf{c}^{T}\mathbf{x}$ decreases yet maintaining all feasibility conditions, hence obtaining an improved basic feasible point.

2.1 Optimality Conditions

At the $k$ -th iteration, let $\left(\mathbf{x}_{k},\mathbf{B}_{k},\mathbf{N}_{k}\right)$ denote the basic feasible point, the constraint coefficient matrix for basic variables $\mathbf{x}_{B}$ and the constraint coefficient matrix for non-basic variable $\mathbf{x}_{N}$ . Let $\mathbf{c}_{B}$ and $\mathbf{c}_{N}$ be cost vector relevant to $\mathbf{x}_{B}$ and $\mathbf{x}_{N}$ . The following conditions need hold in order for $\mathbf{x}_{k}$ to be an optimal solution,

$\displaystyle\mathbf{B}_{k}^{T}\lambda_{B}$	$\displaystyle=$	$\displaystyle\mathbf{c}_{B}$
$\displaystyle\mathbf{N}_{k}^{T}\lambda_{B}+\lambda_{N}$	$\displaystyle=$	$\displaystyle\mathbf{c}_{N}$
$\displaystyle\lambda_{N}$	$\displaystyle\geq$	$\displaystyle 0$

, where $\lambda_{B}$ is the vector of Lagrange multipliers for strict equality conditions $\mathbf{B}_{k}\mathbf{x}_{B}=\mathbf{b}$ , $\lambda_{N}$ is the vector of Lagrange multipliers for non-negative conditions $\mathbf{x}_{N}\geq 0$ where $\mathbf{x}_{N}$ is exactly fixed at 0 currently.

We first solve a $m\times m$ linear system to obtain

\lambda_{B}={\mathbf{B}_{k}^{T}}^{-1}\mathbf{c}_{B}

, and then

\lambda_{N}=\mathbf{c}_{N}-\mathbf{N}_{k}^{T}\lambda_{B}

If all elements in $\lambda_{N}$ are positive, the Simplex algorithm terminates and returns $\mathbf{x}_{k}$ as the solution.

2.2 Pivoting and LU Update

When $\mathbf{x}_{k}$ does not satisfy the optimality conditions, we have $\lambda_{q}<0$ for some $q\in\mathcal{N}$ . Let $\mathbf{A}_{q}$ be the $q$ -th column of matrix $\mathbf{A}$ , which is the constraint coefficient vector for variable $x_{q}$ . It can be shown that there exists a feasible descent direction $\mathbf{r}$ ,

\mathbf{r}=-\mathbf{B}_{k}^{-1}\mathbf{A}_{q}

, such that by moving along this direction,

1.

$x_{q}$ becomes positive and
2.

all other non-basic variables are fixed at 0 and
3.

equality condition $\mathbf{A}\mathbf{x}=\mathbf{b}$ is retained and
4.

objective function value $\mathbf{c}^{T}\mathbf{x}$ decreases.

Also, by moving along $\mathbf{r}$ , if all basic variable $\mathbf{x}_{B}$ increase value, we declare the given LP problem is unbounded. Otherwise, $\exists p\in\mathcal{B}$ such that the basic variable $x_{p}$ hits zero before any other basic varialbes. Now we have an “enter index” $q$ to join the basic index set $\mathcal{B}$ and a “quit index” $p$ to join the non-basic index set $\mathcal{N}$ . This is the pivoting operation. By the end of it, we have obtained a new basic feasible point $\left(\mathbf{x}_{k+1},\mathbf{B}_{k+1},\mathbf{N}_{k+1}\right)$ . One can show that $\mathbf{B}_{k+1}$ obtained by the pivoting operation also have full rank.

A crucial part of a Simplex implementation is to efficiently solve the linear system $\mathbf{B}\mathbf{y}=\mathbf{z}$ for some right hand side vector $\mathbf{z}$ . At the first iteration, we have a LU decomposition on $\mathbf{B}_{0}=\mathbf{L}_{0}\mathbf{U}_{0}$ . Subsequently, $\mathbf{B}_{k+1}$ is obtained by swapping columns $\mathbf{A}_{p}$ and $\mathbf{A}_{q}$ for quit and enter index $p$ and $q$ . The LU decomposition of $\mathbf{B}_{k+1}$ can be economically obtained by updating $\mathbf{L}_{k}$ and $\mathbf{U}_{k}$ . Our implementation is based on the LU update algorithm [4] by J.J. Forrest and J. A. Tomlin.

2.3 Canonical Form Conversion

Given a LP problem descriped in the general form (2), we need to convert it to the canonical form (1). First of all, for an unknown variable $x_{i}$ with simple bound constraint $l_{i}\leq x_{i}\leq u_{i}$ , we can convert it to two equality constraints $x_{i}-\bar{s}_{i}=l_{i}$ and $x_{i}+\hat{s}_{i}=u_{i}$ by introducing slack and surplus variable $\bar{s}_{i}\geq 0$ and $\hat{s}_{i}\geq 0$ . Then $x_{i}$ becomes a free variable. All free variables $x_{i}$ can be further decomposed to positive part $x_{i}^{+}\geq 0$ and negative part $x_{i}^{-}\geq 0$ as $x_{i}=x_{i}^{+}-x_{i}^{-}$ . This approach effectively makes the augmented unknown space very large.

To implement an efficient algorithm, NPL actually only introduces slack and surplus variables to general inequality constraints. For unknowns with simple bound constraints, we construct the index set $\mathcal{B}$ and $\mathcal{N}$ such that the basic variables $\mathbf{x}_{B}$ stay in their feasible region and non-basic variables $\mathbf{x}_{N}$ are fixed at one of their bounds. When we test the optimality conditions, we need flip the sign of a Lagrange multiplier $\lambda_{i}\in\lambda_{N}$ if the non-basic variable $x_{i}$ is fixed at its upper bound.

All general inequality constraints are splitted to a set of upper bounded constraints $\mathbf{A}_{u}\mathbf{x}\leq\mathbf{b}_{u}$ and a set of lower bounded constraints $\mathbf{A}_{l}\mathbf{x}\geq\mathbf{b}_{l}$ . By introducing slack and surplus variables, $\bar{\mathbf{s}}\geq 0$ and $\hat{\mathbf{s}}\geq 0$ , we construct a new set of linear equality constraints with respect to the augmented unknwns $\left(\mathbf{x},\bar{\mathbf{s}},\hat{\mathbf{s}}\right)$ .

\displaystyle\left(\begin{array}[]{ccc}\hidden@noalign{}\hfil\mathbf{A}_{% \mathcal{E}}&&\cr\hidden@noalign{}\hfil\mathbf{A}_{u}&\mathbf{I}&\cr% \hidden@noalign{}\hfil\mathbf{A}_{l}&&-\mathbf{I}\end{array}\right)=\left(% \begin{array}[]{c}\hidden@noalign{}\hfil\mathbf{x}\cr\hidden@noalign{}\hfil% \bar{\mathbf{s}}\cr\hidden@noalign{}\hfil\hat{\mathbf{s}}\end{array}\right)% \left(\begin{array}[]{c}\hidden@noalign{}\hfil\mathbf{b}_{\mathcal{E}}\cr% \hidden@noalign{}\hfil\mathbf{b}_{u}\cr\hidden@noalign{}\hfil\mathbf{b}_{l}% \end{array}\right)

We shall refer this new set of linear constratins as $\mathbf{A}_{aug}\mathbf{x}_{aug}=\mathbf{b}$ . Note that by this construction, we have more unknowns than the linear equality constraints. The constraint specification is then converted to the Canonical form (1).

2.4 Initial Basic Feasible Point

To obtain an initial basic feasible point $\left(\mathbf{x}_{0},\mathbf{B}_{0},\mathbf{N}_{0}\right)$ , we first enforce feasibility to simple bound constraints. If one unknown variable provided by the initial guess is outside of its simple bounds, we fix it at the nearest boundary. Secondly, we add artificial variables $\mathbf{z}\geq 0$ so that equality strictly holds for all linear constraints.

The initial basic index set containts all artificial variables $\mathbf{z}$ and some of $\bar{\mathbf{s}}$ and $\hat{\mathbf{s}}$ . This gives us a diagonal matrix $\mathbf{B}_{0}$ with diagonal elements being either 1 or -1. Lastly, we solve the following phase-I LP feasibility subproblem.

	$\displaystyle\min_{\mathbf{x},\bar{\mathbf{s}},\hat{\mathbf{s}},\mathbf{z}}% \mathbf{e}^{T}\mathbf{z},\mbox{ s.t. }$	$\displaystyle\mathbf{A}_{aug}\mathbf{x}_{aug}\pm\mathbf{z}=\mathbf{b}$
		$\displaystyle\mathbf{l}\leq\mathbf{x}\leq\mathbf{u},$
		$\displaystyle\bar{\mathbf{s}}\geq 0,$
		$\displaystyle\hat{\mathbf{s}}\geq 0,$
		$\displaystyle\mathbf{z}\geq 0,$
		$\displaystyle\mathbf{x}_{aug}=\left(\mathbf{x},\bar{\mathbf{s}},\hat{\mathbf{s% }}\right)$

Clearly, when the Simplex method terminates at the end of this phase-I feasibility problem, we either declare that the original LP problem has no feasible solution if any artificial variable in $\mathbf{z}$ is still positive, or have obtained an initial basic feasible point.

2.5 Selection of the Enter Index $q$

When we test optimality conditions, likely there are multiple negative Lagrange multipliers in $\lambda_{N}$ . There are different strategies in practice to choose the “enter index” $q$ from those negative elements. We use the “Steepest Edge” method outlined in [2], section 13.5. This strategy has the advantage of finding a search direction $\mathbf{r}$ that tends to render large step moves.

3 Active-set Method

Active-set based methods can directly work on the general form (2) of a LP problem. Similar to the Simplex method, an active-set method also generates a sequence of feasible points. At each iteration, it partitions all linear constraints into two sets. One set is called “working set”, which includes the set $\mathcal{E}$ of equality constraints and some inequality constraints that now have equality exactly holds for the current feasible point $\mathbf{x}$ . The other set is the non-working set. Both working and non-working sets are feasible at $\mathbf{x}$ .

Let $\mathcal{W}$ and $\mathcal{N}$ denote the index set of working and non-working general constraints respectively. The size of $\mathcal{W}$ never exceeds $n$ . Let $\mathbf{A}_{\mathcal{W}}$ and $\mathbf{A}_{\mathcal{N}}$ be the corresponding working and non-working constraint coefficient matrix. At each iteration, the method computes a null space direction $\mathbf{p}$ wrt $\mathbf{A}_{\mathcal{W}}$ such that moving along this direction does not change the value of $\mathbf{A}_{\mathcal{W}}$ . Hence we shall refer $\mathbf{p}$ as a feasible direction.

By moving along a descent feasible direction $\mathbf{p}$ , the objective function value will decrease. Non-working constraints will also change values, unless they are linearly dependent to $\mathbf{A}_{\mathcal{W}}$ . If contraint $\mathbf{A}_{i},i\in\mathcal{N}$ hits its bound along $\mathbf{p}$ before any other non-working constraints, strict equality holds for constraint $\mathbf{A}_{i}$ at its bound. It then becomes a working constraint and moves from $\mathcal{N}$ to $\mathcal{W}$ . Whenever the working set $\mathcal{W}$ has exactly $n$ number of constraints, a vertex point is reached, which is a stationary point.

On the other hand, if we can move infinite step length along $\mathbf{p}$ without hitting any non-working constraint, we can declare the LP problem has unbounded solution. This could happen when the rows of $\mathbf{A}_{\mathcal{N}}$ are just linear combinations of $\mathbf{A}_{\mathcal{W}}$ . In other words, there is redundancy in constraint specification.

At a vertex point, optimality conditions are tested on all inequality working constraints $\{i,i\in\mathcal{W}-\mathcal{E}\}$ . If the Lagrange multiplier condition for working constraint $\mathbf{A}_{i}$ is not satisfied, we can move along a negative Lagrange direction such that $\mathbf{A}_{i}$ leaves from its bound and becomes strictly feasible. We also move index $i$ from $\mathcal{W}$ to $\mathcal{N}$ . It can be shown that such a direction is a descent direction. At each iteration, the algorithm repeats to udpate $\mathcal{W}$ and move along a descent feasible direction $\mathbf{p}$ until all optimality conditions are satisfied.

3.1 Simple Bound Constraints

The overall strategy of an active set method is to maintain a set of working constraints and search in a direction that retains the values of these working constraints. The same principle applies to simple bound constraints as well. Without loss of generality, assume all simple bounds are lower bounded constraints $\mathbf{x}\geq\mathbf{l}$ . At each iteration, if a variable is equal to its lower bound, it is considered as a fixed varialbe. A variable that stays stictly inside its feasible region is considered as a free variable.

Let $\mathbf{x}_{b}$ and $\mathbf{x}_{f}$ denote the set of fixed and free variables respectively. We could treat the simple bounds for $\mathbf{x}_{b}$ as general inequality constraints $\mathbf{I}\mathbf{x}_{b}\geq\mathbf{l}_{b}$ with equality exactly being held currently. Therefore we can include as many fixed variables as possible in the current working set and arrange the working constraints in the following block structure,

\left(\begin{array}[]{cc}\mathbf{I}&\cr\mathbf{B}_{\mathcal{W}}&\mathbf{F}_{% \mathcal{W}}\end{array}\right)\left(\begin{array}[]{c}\mathbf{x}_{b}\\ \mathbf{x}_{f}\end{array}\right)=\left(\begin{array}[]{c}\mathbf{l}_{b}\\ \mathbf{b}_{\mathcal{W}}\end{array}\right)

, where the coefficient matrix $\mathbf{A}_{\mathcal{W}}$ of general working constraints is partitioned to $\mathbf{B}_{\mathcal{W}}$ and $\mathbf{F}_{\mathcal{W}}$ for fixed and free variables accordingly. The search direction is a change to those free variables. We just need bookkeep when a free variable hits its bounds and becomes a fixed variable, and when a fixed variable is removed from the working set and becomes free.

3.2 Null Space Direction

Assume there are $m$ number of general constraints in a working set $\mathcal{W}$ . Let $\mathbf{F}_{\mathcal{W}}^{T}$ have a full QR factorization

\mathbf{F}_{\mathcal{W}}^{T}=\underbrace{\left(\begin{array}[]{c|c}\mathbf{Y}&% \mathbf{Z}\end{array}\right)}_{\mathbf{Q}}\left(\begin{array}[]{c}\mathbf{R}% \cr 0\end{array}\right)

The first $m$ columns of $\mathbf{Q}$ forms the range space $\mathbf{Y}$ of $\mathbf{F}_{\mathcal{W}}$ . The remaining columns of $\mathbf{Q}$ is the null space $\mathbf{Z}$ , $\mathbf{F}_{\mathcal{W}}\mathbf{Z}=0$ . If we seek a null space direction of free variables $\mathbf{p}_{f}=\mathbf{Z}\mathbf{p}_{z}$ for some $\mathbf{p}_{z}$ , then $\mathbf{F}_{\mathcal{W}}\mathbf{p}_{f}=0$ . Moving along such a null space direction retains the feasibility with respect to the current working constraints.

To obtain a descent feasible direction, we simply choose $\mathbf{p}_{z}$ being the projection of the negative gradient to the null space, $\mathbf{p}_{z}=-\mathbf{Z}^{T}\mathbf{g}(\mathbf{x}_{f})$ . For LP problems, $\mathbf{g}(x)$ is the cost vector $\mathbf{c}$ . So $\mathbf{p}_{f}=-\mathbf{Z}\mathbf{Z}^{T}\mathbf{c}_{f}$ , where $\mathbf{c}_{f}$ is the cost of free variables. Let $\mathbf{c}_{b}$ denote the cost of fixed variables on the other hand.

3.3 Optimality Conditions

Including the fixed variables as working constraints, the stationary condition requires

\left(\begin{array}[]{cc}\mathbf{I}&\mathbf{B}_{\mathcal{W}}^{T}\cr&\mathbf{F}% _{\mathcal{W}}^{T}\end{array}\right)\left(\begin{array}[]{c}\lambda_{b}\\ \lambda_{\mathcal{W}}\end{array}\right)=\left(\begin{array}[]{c}\mathbf{c}_{b}% \\ \mathbf{c}_{f}\end{array}\right)

Accordingly, we can solve

	$\displaystyle\lambda_{\mathcal{W}}$	$\displaystyle=$	$\displaystyle{\mathbf{F}_{\mathcal{W}}^{T}}^{-1}\mathbf{c}_{f}$
	$\displaystyle\lambda_{b}$	$\displaystyle=$	$\displaystyle\mathbf{c}_{b}-\mathbf{B}_{\mathcal{W}}^{T}\lambda_{\mathcal{W}}$

Next we test Lagrange multipliers $\lambda_{i},i\in\mathcal{W}-\mathcal{E}$ and $\lambda_{i},x_{i}\in\mathbf{x}_{b}$ . For upper bounded constraints, the sign of Lagrange multipliers are flipped. If all $\lambda_{i}$ are positive, the optimal solution is found. Otherwise either the general constraint $\mathbf{A}_{i}$ or the fixed variable $\mathbf{x}_{i}$ that has the most negative $\lambda_{i}$ is selected to quit from the current working set or the set of fixed variables. It can be shown that the feasible direction wrt the updated working set is a descent direction.

3.4 QR Update

The active set method needs repeatedly solve the linear system $\mathbf{F}_{\mathcal{W}}^{T}\mathbf{y}=\mathbf{z}$ for some right hand side vector $\mathbf{z}$ . Given the QR factorization of $\mathbf{F}_{\mathcal{W}}^{T}$ , we solve

\mathbf{R}\mathbf{y}=\mathbf{Y}^{T}\mathbf{z}

$\mathbf{F}_{\mathcal{W}}$ changes whenever a working constraint is added or removed from the working set or the set of free variables changes. An efficient QR factorization update scheme is implemented based on [3], section 5.2.4.

3.5 Initial Working Set and Feasible Point

The active set method also has a phase-I feasibility stage to find an initial feasible point and a subsequent optimality stage to solve the actual LP problem. A “crash start” strategy, [3] section 5.8.1, is implemented to seamlessly connect these two stages. By the end of solving the phase-I feasibility problem, we have not only obtained a pair of initial feasible point and initial working set, but also have all internal resources, i.e. QR factorization, ready to immediately start the the phase-II iterations.

Again, without loss of generality, assume all inequality constraints are lower bound constraints. We first construct a crash index set $\mathcal{C}$ to include all equality constraints $\mathcal{E}$ . Linear dependency is tested so that only independent constraints are added to $\mathcal{C}$ . Next, given an initial guess $\hat{\mathbf{x}}$ and a “crash radius” control parameter $r$ provided by a user, we compute the violation amount $\left(\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}\right)^{+}$ for all general constraints, and $\left(\mathbf{l}-\hat{\mathbf{x}}\right)^{+}$ for simple bounds.

A violated constraint is added to $\mathcal{C}$ if its violation amount is smaller than $r$ and it is independent to all existing constraints in $\mathcal{C}$ . This procedure stops as soon as $\mathcal{C}$ has $n$ constraints.

Next, we can fix all violated variables in $\mathcal{C}$ to the corresponding simple bound values and obtain $\mathbf{x}_{b}$ . Now let $\mathcal{W}$ include all general constraints in $\mathcal{C}$ and solve free variables from

\mathbf{F}_{\mathcal{W}}\mathbf{x}_{f}=\mathbf{b}-\mathbf{B}_{\mathcal{W}}% \mathbf{x}_{b}

At this point, we have obtained an initial working set $\mathcal{W}$ and an initial point $\left(\begin{array}[]{c}\mathbf{x}_{b}\cr\mathbf{x}_{f}\end{array}\right)$ that is feasible to working general constraints $\mathbf{A}_{\mathcal{W}}$ . Of course, some of the non-working general constraints $\mathbf{A}_{\mathcal{N}}$ or simple bounds could still be violated. Conceptually, we introduce artificial variables $\mathbf{y}$ and $\mathbf{z}$ and solve a phase-I feasibility subproblem,

	$\displaystyle\min_{\mathbf{x},\mathbf{y},\mathbf{z}}\mathbf{e}^{T}\mathbf{y}+% \mathbf{e}^{T}\mathbf{z},\mbox{ s.t. }$	$\displaystyle\mathbf{b}_{l}\leq\mathbf{A}\mathbf{x}\pm\mathbf{z}\leq\mathbf{b}% _{u}$
		$\displaystyle\mathbf{l}\leq\mathbf{x}\pm\mathbf{y}\leq\mathbf{u},$
		$\displaystyle\mathbf{y}\geq 0$
		$\displaystyle\mathbf{z}\geq 0$

When the active set method terminates after its phase-I stage, if the objective function value is not zero, we declare the LP problem has no feasible solution. Otherwise, the solution and working set after phase-I is the starting point of phase-II.

4 Termination Control

Both the Simplex method and the active set method have a control parameter $N$ for the max number of iterations allowed. The algorithm always terminates after $N$ iterations. Further, there is an exception-throw control. Once enabled, the algorithm throws an exception after $N$ iterations. If disabled, it just reports the “half-way” solution found after $N$ iterations.

Both algorithms also have a “convergence tolerance” parameter $\epsilon_{c}>0$ that is used to check the Lagrange multipliers in optimality condition test. If the most negative Lagrange multiplier $\lambda\geq-\epsilon_{c}$ , the Lagrange multiplier condition is considered as satisfied.

In Simplex method, the number of basic variables is exactly equal to the number of equality constraints. So the basic feasible points generated by the algorithm are always vertex points. Stationary condition is automatically satisfied at a vertex point.

In the active set method, the number of working constraints usually is equal to the number of variables $n$ , hence a vertex point as well. However, it is possible that the cost vector $\mathbf{c}_{f}$ for free variables happens to be 0 or a linear combination of $\mathbf{F}_{\mathcal{W}}$ at a non-vertex point. The active set method has a “stationary tolerance” control parameter $\epsilon_{s}>0$ . When $|\mathbf{c}_{f}|_{\infty}\leq\epsilon_{s}$ or $|\mathbf{c}_{f}^{T}\mathbf{p}_{f}|\leq\epsilon_{s}$ , we consider the stationary condition is satisfied at a non-vertex point.

5 Summary

The active set method is a very powerful tool to handle LP and many other contrained optimization problems. It directly takes the general LP form (2) without having to introduce slack and artifical variables. When the number of general contraints $m$ is significantly greater than the number of unknowns $n$ , it has the advantage of working on a linear system that is no bigger than the number of free variables.

Simplex method simultaneously works on all linear constraints. Initially, we need convert the general form to canonical form (1) with slack and artifical variables. It is most suitable for the situation $n\gg m$ .

References

[1] David G. Luenberger: Linear and Nonlinear Programming, 2nd edition, Springer 2003
[2] Jorge Nocedal and Stephen J. Wright: Numerical Optimization, Springer, 1999
[3] Philip E. Gill, Walter Murray and Margaret H. Wright: Practical Optimization, Academic Press, 1981
[4] J.J. Forrest and J. A. Tomlin: Updated Triangular Factors of the Basis to Maintain Sparsity in the Product from Simplex Method, Mathematical Programming 2(1972) 263-278

Linear Programming Mathwrist White Paper Series