LP  Simplex Draft
Contents
Everybody solves some optimization problems – the airlines schedule flights, companies manage production facilities, salesman still looks for traveling options… When you need to achieve the best outcome minimizing/maximizing a linear cost functions you meet linear programming.
LP was developed for military purposes in 1939, so it has a long history. Now it is one of the most important problems in operations research and heavily used in different areas directly or as subproblems.
Intro
The standard algorithm here is the simplex method. Lots of information is floating around, but let’s review the basics. There’re several forms for expressing linear programs. Standard form:
$$\begin{array}{1}
\text{maximize }{c}^\prime x \\
\text{subject to}\\
Ax <= b\\
x >= 0
\end{array}$$
We assume that the rows of matrix $A$(m x n) are linearly independent.
The form above can be easily converted to slack form with a help of the new variables (slack variables):
$$\begin{array}{1}
s = b_i  \sum a_{ij} x_j \\
s >= 0
\end{array}$$
So the problem now is
$$\begin{array}{1}
\text{maximize }{c}^\prime x\\
\text{subject to}\\
Ax = b\\
x >= 0
\end{array}$$
Geometrical interpretation for 2 variables is very intuitive: all $x_1$ and $x_2$ satisfying the constraints are feasible solutions, and an optimal solution – maximum or minimum – is at a vertex. This feasible region, for n variables – in ndimensional space, is called a simplex. The algorithm terminates when it reaches an optimal objective value at some vertex.
Why write a custom implementation? For me Excel solver is usually a solution ^{1}, works like a charm – when it’s not a Mac version. The naïve version of algorithm is quite simple to implement, it’s also a nice refresher and managed to solve the problem, so why not?
Consider the following problem:
$$\begin{array}{1}
\text{minimize }x_1 + 5x_2  2x_3 \\
\text{subject to} \\
x_1 + x_2 + x_3 <= 4 \\
x_1 <= 2 \\
x_3 <= 3 \\
3x_2 + x_3 <= 6 \\
x_1, x_2, x_3 >= 0
\end{array}$$
The solution is very straightforward:
1. Convert to slack form
$$\begin{array}{l}
\text{minimize } x_1 + 5x_2  2x_3\\
\text{subject to }\\
x_1 + x_2 + x_3 + x_4 <= 4\\
x_1 + x_5 <= 2\\
x_3 + x_6 <= 3\\
3x_2 + x_3 + x_7 <= 6\\
x_1,\ x_2,\ x_3,\ x_4,\ x_5,\ x_6,\ x_7 <= 0
\end{array}$$
2. We see that the cost can be reduced with increasing $x_3$ (it’s also obvious that max $x_3$ value is equal to 3). If the cost is already optimal the solution is found.
For each $j$ the reduced cost is defined as
$$\begin{array}{1}
\bar{c_j} = c_j  {c_b}^\prime B^{1}A_j\\
\text{where }c_b\text{ – the vector of basic variables costs.}
\end{array}$$
3. Look for a direction of cost decrease, when moving into the new direction we replace a basic variable with a new one; go to the next iteration.
Implementation
The first version I came up with was entirely immutable – with all the costs of creating new sets, matrices etc. But there’s no sense to recreate the whole matrix when only one column is changed for the next iteration. This version is listed below – not that functional, more dangerous, with a mutable state, but cares about time/memory. Note, that this implementation doesn’t handle some cornercases as it was not a part of my goal.


For simplicity we use the naïve initialization – all given variables become nonbasic, the slack ones – basic and their values are equal to constraints vector $$b$$. The full method returns the solution vector and cost function value:


And now time for my favorite part of the method and one the best ways to improve performance – reuse the data you already have to avoid recomputations whenever it’s possible. In this case it’s about inverting the basis matrix. We know that $B$ at the next iteration is the same as current, except one column. There’s also the current $B^{1}$. The naïve simplex + math = revised simplex method.


Sample


Complete snippet version is available at github.
Summing up
 mutability may be evil, but sometimes it’s a way to go;
 reusing the computation results and the parts of datastructures ftw.
What to look at
 MS Solver Foundation – quite a nice tool to check too (with msi downloads, but I tried it a long time ago with mono – and it worked).
 some theory: Introduction to Linear Optimization book (Dimitris Bertsimas and John N. Tsitsiklis) or anything else.
 one of the sideeffects of my work are numerous spreadsheets, and Solver proved to be incredibly helpful. ^{[return]}
Author Natallie
LastMod 20130622