Simplex

Simplex is the first systematically method for solving linear program via a sequence of Gaussian eliminations. It was invented by George Dantzig in 1947 [1]. Algebraically, it converts the linear program from one slack form into to an equivalent one whose objective value does not decrease and possibly increase. It keeps going until the optimal is reached. Geometrically, it is a greedy search strategy that moves from a vertex of the feasible region to another one, with the objective value at the new vertex as least as good as or possibly improved over the previous one. It is an exponential algorithm, although most of time it does a lot better. In this blog we walk through a simple example to illustrate this.

In this case we are operating a factory with two products. Both product have 1 unit of profit. Making the products needs to use two kind of machine. Product one uses 1 hour of machine one and product two uses 2 hours of machine one. Further, product one uses 2 hours of machine two and product two uses 1 hour of machine two. Both machines operate at most 6 hours each day. We want to maximize our profit:

It can be represented concisely in matrix form:

where

We call the set of points that satisfy the constraints of the optimization problem the feasible region

It is obvious that the optimal is obtained at $(2, 2)$, which gives the objective value $4$. It is a vertex of the feasible region.

Slack Form

In general, it might not be obvious what the optimal point is. But if the elements of $b$ are non-negative, the origin is a feasible point, from which we can begin our search. To make the search easier, we are going to introduce an additional number of variables, the slack variables $s_1, s_2$ and $z$ that correspond to the constraints and objective value.

rewrite the optimization into slack form:

They are called slacks variables because they correspond how much slack you have in the inequalities in the original problem. For this slack form, we also call $s_1$ and $s_2$ basic variables, and the original variables non-basic variables.

Now the feasible region $F$ is considered as the intersection of the two subspaces

We can represent the LP in slack form even more concisely in a table:

where the last column denotes the objective function. This example has a trivial starting point : $(0, 0, 6, 6)$. It is called a basic solution in which all non-basic variables are set to zero. It can be read from the tabular form: the identity matrix correspond to basic variables, and $(s_1, s_2) = b$.

Pivoting

Look at the original LP:

As we would like to maximize the objective, we would like to increase the value of $x_1$ or $x_2$ as much as possible. Suppose we just increase one variable at a time. For example, we can increase variable $x_1$ and keep the value of $x_2$ to be $0$.

How much can we increase $x_1$? After increasing $x_1$, we need to maintain the inequality constraints:

The values of $s_1$ and $s_2$ decrease as $x_1$ increases. By they can be drop below $0$. The maximum possible increases of $x_1$ is bounded by $6 / 2 = 3$. If we set $x_1 = 3$, $s_1 = 6 - x_1 = 3$ and $s_2 = 6 - 2x_1 = 0$.

Note that $x_1$ increases from $0$ to $3$ and $s_2$ decreases from $6$ to $0$. This process is called pivoting.

How do we perform pivoting in the tabular form? By the replacement we have just performed -- normalizing the second row, and using it to eliminate all other elements in the first column.

Recall the original table:

We divide the second row by a factor of 2 to normalize it, and subtract the first and third row by the second row

Observe the first two rows of the tabular form:

If we set $x_2 = 0$ and $s_2 = 0$, the constraints they implies becomes

Again the unity columns correspond to non-zero elements and we can read the solution from the tabular form.

Correctness of Pivoting

Before we continue our search on the tabular form, we need to show that pivoting is correct in the sense that, after pivoting, the tabular form represents the same optimization.

In particular, we need to guarantee that the feasible regions are the same and the objective function are equivalent defined on the feasible regions.

We focus on the extended matrix that represents the constraints:

FACT 1. Doing a row replacement on $\bar A$ does not change $F$, where row replacement is defined as

i.e., adding $k$ times the $j$-th row to the $i$-th row, for $i \neq j, k \in R$.

To see this, denote $F$ and $F'$ the sets before and after the replacement respectively. $F = F'$ implies that each point $x$ belonging to $F$ belongs to $F'$, and vice-versa. Specifically, the row replacement corresponds to left multiplying $\bar A$ by a matrix

which is a diagonal matrix with an additional $k$ in the $(i, j)$ position. It is invertible with inverse matrix

The previous claim becomes $\{ x : E_{(i, j), k} A x = E_{(i, j), k} b \} = \{ x : A x = b \}$.

But wait. Why do we subtract the last row from the second one? Fact 1 has justified the replacement operation on $\bar A$, but not on the objective coefficients $c$ (the last row in the tabular form). To illustrate the meaning of the operation, recall that

Hence $z$ is the value of the objective function. After normalizing the second row, it becomes

The constraint holds for any feasible point of the linear programming. Therefore,

The last $-3$ in the table is the inverse constant in the objective functions after the replacement of variable.

Back to our search

To improve the current solution further, we note that the coefficient of $x_2$ is positive. We can increase the value of $x_2$. The maximum increase is given by $\min \{ 3 / (1 / 2), 3 / (3/ 2) \} = 2$, determined by the first row:

To make the change explicitly, we perform elimination by the first row

which give solution $x = (2, 2, 0, 0)$ and objective value $4$. It is indeed optimal solution for the problem, as all coefficients of the optimization are non-positive.

Optimality Condition and Strong Duality

We investigate the result deeper by viewing the initial tabular form in abbreviation as

The series of $m$ Gaussian elimination can be interpreted as left multiplying a matrix:

where each $E_j$ ( $1 \le j \le m$) matrix is invertible and represents a Gaussian elimination, i.e., if we left multiple a matrix by $E_j$, it is equivalent to perform Gaussian elimination on the matrix. The matrix product $E_m...E_3 E_2 E_1$ has the form

because we never use the last row to eliminate the other rows.

Now:

In our example,

we have

The optimality condition of simplex algorithm implies that

This is the optimal solution for the dual program of the primal:

If simplex algorithm stops, it implies strong duality holds:

• The optimal value of the primal program equals to the optimal one of the dual program, assuming that both values exist and are bounded.

Geometric Interpretation

In the beginning we claim that the simplex algorithm move from one vertex to another. We prove it rigorously in this section. First we need a few definitions:

Vertex

A vertex is a point $v \in F$, such that $\nexists v_1, v_2 \in F, v_1 \neq v_2$ and $\lambda \in (0, 1)$, s.t., $v = \lambda v_1 + (1 - \lambda) v_2$.

In other words, $v$ is a vertex of $F$ if it is not contained in a line segment of $P$.

Theorem. If LP has an optimal solution, then it has an optimal solution that is a vertex of its feasible set.

Proof. Denote $v$ the optimal solution of the LP with maximum number of zero components. If $v$ is not a vertex, then $\exists v_1, v_2 \neq v, v_1, v_2 \in F$ and $\lambda > 0$, such that $\lambda v_1 + (1 - \lambda) v_2 = v$.

As $v$ is optimal, we have $c^T v \ge c^T v_1$ and $c^T v \ge c^T v_2$. By $\lambda v_1 + (1 - \lambda) v_2 = v$, it concludes that $c_T v = c^T v_1 = c^T v_2$.

Define $I = \{ i : v[i] > 0 \}$. As $v_1 \ge 0, v_2 \ge 0$ and $\lambda v_1 + (1 - \lambda) v_2 = v$ so that for $i \neq I$, $v_1[i] = v_2[i] = v[i] = 0$.

But $A_I v_1[I] = A_I v_2[I] = b$. Let $x = \epsilon (v_1 - v_2) + v$.

1. $x$ is feasible, for small enough $\epsilon$.
2. $x$ is optimal, since $c^T x= \epsilon c^T (v_1 - v_2) + c^T v = c^T v$.
3. $x_i = 0$ for $i \notin I$.

We can find some $\epsilon$ (either positive or false) to make $x_i = 0$ for some $i \in I$. Now $x$ contains more zeros than $v$, a contradiction. $\square$

Basic Solution

A solution $x$ of $Ax = b$ is called a basic solution if $\{A_i : x_i \neq 0 \}$ (columns in $A$ that correspond to non-zero components in $x$ ) are linearly independent.

Basic Feasible Solution

A basic solution $x$ is called basic feasible solution if $x \ge 0$, i.e., $x \in \{ x : Ax = b\} \cap \{ x \ge 0 \}$.

With respect to the feasible region of LP in slack form, the vertices have have following property:

Lemma:. A point $v \in F = \{x : Ax = b\} \cap \{ x: x \ge 0\}$ is a vertex of $F$ if and only if it is a basic feasible solution. Here $A \in R^{m \times n}, x \in R^n, b \in R^{m}$ and $rank(A) = m$ (otherwise some row constraints of the $Ax = b$ are redundant).

Proof: Denote $S$ the set of indices of positive components in $v$ and $A_S$ the set of

Only if: If $A_S$ does not have full column rank, then $\exists u \in R^{|S|}$, such that $A_S u = 0$. Then for small enough $\epsilon \in R, \epsilon > 0$, both $v + \epsilon u$ and $v - \epsilon u$ are feasible point in $F$. Now $v = 0.5(v + \epsilon u) + 0.5(v - \epsilon u)$, a contradiction.

If: As columns of $A_S$ are linearly independent, the only solution of to $A_S x_S = 0$ is given by $x_S = 0$. Hence, $v_S$ is now the unique solution to $A_S v_S = b$. If $v$ is not a vertex, $\exists v_1, v_2, \lambda$, such that $v = \lambda v_1 + (1 - \lambda) v_2$. But in this case we must have $v_1 = v_2 = v$, as $(v_1)_S = (v_2)_S = v_S$ and $(v_1)_{\bar S} = (v_2)_{\bar S} = v_{\bar S} = 0$. A contradiction.

$\blacksquare$.

Corollary: The simplex algorithm moves from vertex to vertex.

Proof: The positive components in the solution correspond to the columns in the identity matrix in the tabular form of simplex, which are linearly independent.

Correctness of Simplex (TO FINISH)

Though not shown in the example, there are two issue associated with simplex algorithm.

1. How do we find an initial solution? If there are negative component in $b$, then $x = 0 \wedge s = b$ is not a feasible solution.
2. How do we guarantee that simplex terminates?

In the previous example, at each step we choose to increase the variable with largest coefficient. In some rare cases, this results in a loop of the algorithm.

One method for avoiding looping is called Bland's rule:

1. If there is a positive coefficient in the objective function, choose the one with largest coefficient.
2. If there are multiple rows with tight constraints, choose the one with largest coefficient.

Theorem. Blands' rule guarantees that the algorithm stops.

Proof. Suppose that on the contrary, the algorithm loops. Denote $B_1, B_2, ..., B_k$ the set of basis in the loop. In each base, there is entering variable and one leaving variable. Every variable that leaves the base must enter another base later. We called these variable fickle variables. Denote $x_t$ the fickle variable with the largest index.

1. $D$ the dictionary that $x_t$ leaves the base, and $x_e$ the variable that enters the base.
2. $D'$ the dictionary that $x_t$ enters the base again.

We have

1. $x_e$ is the entering variable in $D$, therefore $c_e > 0$.
2. $x_t$ is the entering variable in $D'$, and $e < t$, by the definition of $x_t$, we have $x_e' \le 0$.
3. Hence, $y^T A_N[:, e] = c_e' - c_e < 0$.
4. $\exists i, s.t., y_i A_{i, e} < 0$. Then $y_i \neq 0$, and $x_i$ is fickle. By definition, $i < t$. But $x_i$ is not entering $D'$, we know that $y_i < 0$.
5. $A_{i, e} > 0$. But since $x_i$ is fickle, we have $b_i = 0$. In this case, $x_i$ is chosen in preference to $x_t$. A contradiction.

Reference

[1]. Dantzig, George B. "Origins of the simplex method." A history of scientific computing. 1990. 141-151.