This blog is talking about the naïve implementation of simplex method and then focus on how to speed up simplex method. And the second idea induce the Simplex Tableau implementation.

Now, we analysis the computation complexity in current simplex algorithm.

1. Inverting matrix $[B]_{m \times m}$ takes $O(m^3)$ arithmetic operations.

2. Solving an $m \times m$ linear system takes $O(m^3)$ arithmetic operations.

3. Computing $[B][\bar{b}]$ takes $O(m^2)$ arithmetic operations.

4. Computing $\bar{p}^\top \bar{b}$ takes $O(m)$ arithmetic operations.

Naïve Implementation

At the begin of iteration, we have the indices $B(1), B(2), …, B(m)$ of the current basic variables and the basis matrix $B$.

We compute $\bar{p}^\top := c_B^\top B^{-1}$, which consumes $O(m^3)$ operations. The vector $\bar{p}$ is called the vector of simplex multipliers associated with the basis matrix $B$.

Then the reduced cost $\tilde{c}_i = c_i - c_B^T B^{-1} A_i = c_i - \bar{p}^\top A_i$, which takes $O(mn)$ operations.

Find a (non-basic) variable $x_j$ with negative reduced cost $\tilde{c}_j < 0$. The corresponding column $A_j$ is selected to enter the basis.

Solve $B\bar{\mu} = \bar{A_j}$ to find $\bar{\mu} = B^{-1} \bar{A}_j$ from which we get the direction $\bar{d}$ along which we will be moving away from the current basis feasible solution. (But, here we don’t need to compute $B^{-1}$ again).

\[\bar{d} = \begin{bmatrix} \bar{d}_B \\ \bar{d}_N \end{bmatrix} = \begin{bmatrix} - \bar{\mu} \\ \bar{e}_j \end{bmatrix}\]

Finally, determine $\theta^*$ and the variable $x_{B(\ell)}$ that will exit the basis and construct the new basis solution.

\[\theta^* = \min\limits_{i \in \{1,2,\cdots,m\} \text{ s.t. } \mu_i > 0} \left( \frac{x_{B(i)}}{u_i} \right) = \frac{x_{B(\ell)}}{\mu_{\ell}}\]

Total computational effort of a simplex iteration: $O(m^3 + mn)$

Revised Implementation

In this revised implementation we can finish each loop only in $O(m^2 + mn)$.

At begin of iteration, make the matrix $B^{-1}$ available. Then compute $\bar{p}^T = \bar{c}_B^\top B^{-1}$ and $\bar{\mu} = B^{-1} \bar{A}_j$. Now will find another way updating the matrix $B^{-1}$ each time more efficient.

Basis matrix: \(B = \begin{bmatrix} \bar{A}_{B(1)} & \cdots & \bar{A}_{B(\ell-1)} & \bar{A}_{B(\ell)} & \bar{A}_{B(\ell+1)} & \cdots & \bar{A}_{B(m)} \end{bmatrix}\)

Updated basis matrix: \(B' = \begin{bmatrix} \bar{A}_{B(1)} & \cdots & \bar{A}_{B(\ell-1)} & \bar{A}_j & \bar{A}_{B(\ell+1)} & \cdots & \bar{A}_{B(m)} \end{bmatrix}\)

Key idea: Use knowledge of $B^{-1}$ to simplify the computation of $B’^{-1}$.

Definition: Any operation on a matrix $M$ of type “replace the $i^{th}$ row of $M$ by itself plus $\beta$ times the $j^{th}$ row of $M$” is called an elementary row operation on the matrix $M$.

We see that multiplying the $j^{th}$ row by $\beta$ and adding it to the $i^{th}$ row (for $i \neq j$) is the same as left-multiplying by the matrix $Q = I + D_{ij}$, where $D_{ij}$ is a matrix with all entries equal to zero, except for the $(i,j)\text{th}$ entry which is equal to $\beta$. The determinant of such a matrix $Q$ is equal to 1 and, therefore, $Q$ is invertible.

Suppose now that we apply a sequence of $K$ elementary row operations and that the $k^{th}$ such operation corresponds to left-multiplication by a certain invertible matrix $Q_k$. Then, the sequence of these elementary row operations is the same as left-multiplication by the invertible matrix $Q_KQ_{K-1}\cdots Q_2Q_1$. We conclude that performing a sequence of elementary row operations on a given matrix is equivalent to left-multiplying that matrix by a certain invertible matrix.

Since $B^{-1}B = I$, we see that $B^{-1}\bar{A}_{B(i)}$ is the $i^{th}$ unit vector $\bar{e}_i$. Using this observation, we have

\[B^{-1}B = \begin{bmatrix} | & | & | & | & | \\ \bar{e}_1 & \cdots & \bar{e}_{\ell-1} & \bar{\mu} & \bar{e}_{\ell+1} & \cdots & \bar{e}_m \\ | & | & | & | & | \end{bmatrix} = \begin{bmatrix} 1 & & \mu_1 & & \\ & \ddots & \vdots & & & \\ & & \mu_{\ell} & & \\ & & \vdots & \ddots & &\\ & & \mu_{\ell} & & 1 \end{bmatrix}\]

where $\bar{\mu} = B^{-1}\bar{A}_j$. Let us apply a sequence of elementary row operations that will change the above matrix to the identity matrix. In particular, consider the following sequence of elementary row operations.

(a) For each $i \neq \ell$, we add the $\ell^{th}$ row times $-\mu_i/\mu_{\ell}$ to the $i^{th}$ row. (Recall that $\mu_{\ell} > 0$.) This replaces $u_i$ by zero.

(b) We divide the $\ell^{th}$ row by $\mu_{\ell}$. This replaces $\mu_{\ell}$ by one.

Performing this sequence of elementary row operations is equivalent to left-multiplying $B^{-1}B’$ by a certain invertible matrix $Q$ so as to obtain $I$:

\[Q \cdot B^{-1}B' = I\] \[Q \cdot B^{-1} = \left(B'\right)^{-1}\]

Conclusion: Performing the elementary row operations (a) and (b) on $B^{-1}$ returns $\left(B’\right)^{-1}$.

Iteration of the revised simplex algorithm

1. In a typical iteration, we start with a basis consisting of the basic columns $\bar{A}{B(1)}, \cdots , \bar{A}{B(m)}$, an associated basic feasible solution $x$, and the inverse $B^{-1}$ of the basis matrix.

2. Compute the row vector $p^\top = c_B^\top B^{-1}$ and then compute the reduced costs $\tilde{c}_j = c_j - p^\top \bar{A}_j$. If they are all nonnegative, the current basic feasible solution is optimal, and the algorithm terminates; else, choose some $j$ for which $\tilde{c}_j < 0$.

3. Compute $\mu = B^{-1} \bar{A}_j$. If no component of $\mu$ is positive, the optimal cost is $-\infty$, and the algorithm terminates.

4. If some component of $\mu$ is positive, let

\[\theta^* = \min_{\{i=1,...,m \mid \mu_i>0\}} \frac{x_{B(i)}}{\mu_i}\]

1. Let $\ell$ be such that $\theta^* = x_{B(\ell)}/\mu_{\ell}$. Form a new basis by replacing $\bar{A}_{B(\ell)}$ with $\bar{A}_j$. If $y$ is the new basic feasible solution, the values of the new basic variables are $y_j = \theta^*$ and $y_{B(i)} = x_{B(i)} - \theta^*\mu_i$, $i \neq \ell$.

2. Form the $m \times (m + 1)$ matrix $[B^{-1} \mid \bar{\mu}]$. Add to each one of its rows a multiple of the $\ell^{th}$ row to make the last column equal to the unit vector $\bar{e}_{\ell}$. The first $m$ columns of the result is the matrix $(B’)^{-1}$.

Tableau Implementation

Now, we will give an implementation of simplex method in tableau. Here, instead of maintaining and updating the matrix $B^{-1}$, we maintain and update the $m \times (n + 1)$ matrix

\[B^{-1}\begin{bmatrix} \bar{b} \mid A \end{bmatrix}\]

with columns $B^{-1}\bar{b}$ and $B^{-1}\bar{A}_1, …, B^{-1}\bar{A}_n$. This matrix is called the simplex tableau. Note that the column $B^{-1} \bar{b}$, called the zeroth column, contains the values of the basic variables. The column $B^{-1} \bar{A}_i$ is called the $i^{th}$ column of the tableau. The column $\mu = B^{-1} \bar{A}_j$, corresponding to the variable that enters the basis is called the pivot column. If the $\ell^{th}$ basic variable exits the basis, the $\ell^{th}$ row of the tableau is called the pivot row. Finally, the element belonging to both the pivot row and the pivot column is called the pivot element. Note that the pivot element is $\mu_{\ell}$ and is always positive (unless $\mu \leq 0$, in which case the algorithm has met the termination condition in Step 3).

Now, we will introduce how to update the tableau $B^{-1}\begin{bmatrix} \bar{b} \mid A \end{bmatrix}$. This can be accomplished by left-multiplying the current tableau by the matrix of the revised simplex tableau with a matrix $Q$ satisfying $QB^{-1} = (B’)^{-1}$.

Recall the revised simplex, $x_{B(i)} / \mu_i$ is the ratio of the entry in the zeroth column of the tableau to the $i^{th}$ entry in the pivot column of the tableau. We only consider those $i$ for which $u_i$ is positive. The smallest ratio is equal to $\theta^*$ and determines pivot row $\ell$.

Then we give the structure of the tableau as below:

The entry at the top left corner contains the value $-c_B^\top x_B$, which is the negative of the current cost. (The reason for the minus sign will be explained in following) The rest of the zeroth row is the row vector of reduced costs, that is, the vector $\tilde{c} = \bar{c} - c_B^\top B^{-1}A$.

$ -c_B^\top B^{-1}\bar{b} $	$ \bar{c} - c_B^\top B^{-1}A $
$ B^{-1} \bar{b} $	$ B^{-1}A $

Updating rule for zeroth row

First, let’s look at the top left corner of the full tableau, i.e., the opposite of the cost of the current b.f.s.

\[\begin{aligned} - \bar{c}_B^\top x'_B = - c^\top \bar{x}' &= - \bar{c}^\top (\bar{x} + \theta^* \bar{d})\\ &= - (\bar{c}^\top \bar{x} + \theta^* c^\top \bar{d} )\\ &= - (\bar{c}_B^\top \bar{x}_B + \frac{x_{B(\ell)}}{\mu_{\ell}} \tilde{c}_j ) \\ & = - \bar{c}_B^\top \bar{x}_B - \frac{\tilde{c}_j}{\mu_{\ell}} x_{B(\ell)} \end{aligned}\]

Next, let’s look at the change in the simplex multipliers row vector:

\[\begin{aligned} \bar{c}_{B'}^\top (B')^{-1} - \bar{c}_{B'}^\top (B^{-1}) &= \bar{c}_{B'}^\top Q B^{-1} - \bar{c}_{B}^\top B^{-1} = (\bar{c}_{B'}^\top Q - \bar{c}_{B}^\top) B^{-1}\\ &= (\begin{bmatrix} c_{B'(1)} & \cdots & c_{B'(\ell)} & \cdots & c_{B'(m)} \end{bmatrix} Q - \bar{c}_B^\top ) B^{-1}\\ &= (\begin{bmatrix} c_{B(1)} & \cdots & c_{B(\ell -1)} & \blacksquare & c_{B(\ell + 1)} & \cdots & c_{B(m)} \end{bmatrix} - \bar{c}_B^\top ) B^{-1}\\ &= \begin{bmatrix} 0 & \cdots & 0 & * & 0 & \cdots & 0 \end{bmatrix} B^{-1} \end{aligned}\]

where, we know

\[\begin{aligned} * &= - \frac{1}{\mu_{\ell}} (c_{B(1)} \mu_1 + \cdots + c_{B(\ell-1)} \mu_{\ell-1} - c_j + c_{B(\ell+1)} \mu_{\ell+1} + \cdots + c_{B(m)} \mu_m) - c_{B_{(\ell)}}\\ & = \frac{c_j}{\mu_{\ell}} - \frac{1}{\mu_{\ell}} (c_{B(1)} \mu_1 + \cdots - c_{B_{(\ell)}} \mu_{\ell} + \cdots + c_{B(m)} \mu_m) \\ & = \frac{1}{\mu_{\ell}} ( c_j - \bar{c}_B^\top \bar{\mu} ) \\ &= \frac{1}{\mu_{\ell}} (c_j - \bar{c}_B^\top B^{-1} \bar{A}_j) \\ & = \frac{-\tilde{c}_j}{\mu_{\ell}} \end{aligned}\]

It follows that

\[\bar{c}^\top - \bar{c}^\top_B (B')^{-1} A = \bar{c}^\top - \bar{c}^\top_B B^{-1} A - \frac{\tilde{c}_j}{\mu_{\ell}} \bar{e}_{\ell}^\top B^{-1} A\]

Conclusion: The rule for updating the zeroth row turns out to be identical to the rule used for the other rows of the tableau: add a multiple of the pivot row to the zeroth row to set the reduced cost of the entering variable to zero.

Iteration of the full tableau implementation

1. A typical iteration starts with the tableau associated with a basis matrix $B$ and the corresponding basic feasible solution $x$.

2. Examine the reduced costs in the zeroth row of the tableau. If they are all nonnegative, the current basic feasible solution is optimal, and the algorithm terminates; else, choose some $j$ for which $\tilde{c}_j < 0$.

3. Consider the vector $\mu = B^{-1} \bar{A}_j$, which is the $j^{th}$ column (the pivot column) of the tableau. If no component of $\mu$ is positive, the optimal cost is $-\infty$, and the algorithm terminates.

4. For each $i$ for which $u_i$ is positive, compute the ratio $x_{B(i)}/\mu_i$. Let $\ell$ be the index of a row that corresponds to the smallest ratio. The column $\bar{A}_{B(\ell)}$ exits the basis and the column $\bar{A}_j$ enters the basis.

5. Add to each row of the tableau a constant multiple of the $\ell^{th}$ row (the pivot row) so that $\mu_{\ell}$ (the pivot element) becomes one and all other entries of the pivot column become zero.