# Solving LPs graphically and by brute-force using Python

In order to understand better the properties of Linear Programs (LP), it can be helpful to look at some naive methods. In this post, I will solve the following LP graphically using Matplotlib and then by brute-force using NumPy.

## Graphically

Let us first rewrite the constraints.

We have now three functions we can plot.

To find the feasible region we can look at the generated plot. The black lines were drawn manually using GIMP. The black region defines our convex polyhedron. The fundamental theorem of linear programming tells us now that the corners are the maxima of the region. This is also intuitive, because these are the only points that cannot be written as a convex combination of other points (convex independent).

The extreme points are then $E = \left\{\begin{bmatrix}0\\0\end{bmatrix}, \begin{bmatrix}0\\4\end{bmatrix}, \begin{bmatrix}5\\0\end{bmatrix}, \begin{bmatrix}4\\2\end{bmatrix}, \begin{bmatrix}1\\5\end{bmatrix}\right\}$.

If we plug the vectors into the objective function $f\left(\mathbf{x}\right) = x_1 + x_2$, we will find that $f\left(\begin{bmatrix}4\\2\end{bmatrix}\right) = f\left(\begin{bmatrix}1\\5\end{bmatrix}\right) = 6$.

The maximum of $f$ defined over our polyhedron is thus $6$. This value occurs at the points $\text{Conv}\left(\left\{\begin{bmatrix}1\\5\end{bmatrix}, \begin{bmatrix}4\\2\end{bmatrix}\right\}\right)$ where $\text{Conv}$ defines the convex hull of a set.

## Brute-force

The second way to solve an LP uses the algebraic definition of extreme points. This time the corner points are obtained by solving a linear system of $m$ equations. Let’s rewrite the LP in matrix form.

Since the first three constraints are in $\mathbb{R}^2$, there are $\binom{m}{n} = \binom{3}{2} = 3$ possible extreme points / intersections. A simple algorithm would be then

1. Choose $n$ constraints (rows) from the $m \times n$ matrix $\mathbf{A}$.
2. Check if the constraints are linearly independent (e.g. $x - y = 2$ and $x - y = 0$ would be parallel).
3. Solve the system of equations for $\mathbf{x}$.
4. If $\mathbf{x}$ satisfies all constraints, calculate $\mathbf{c}^T\mathbf{x}$, save result and go back to 1.

It is often easier to deal with LPs in standard form, because the matrix inequality gets replaced by $\mathbf{A}\mathbf{x} = \mathbf{b}$. The matrix becomes $% $. Before we had an overdetermined system (more equations than unknowns), now $\mathbf{A}$ is underdetermined (more unknowns than equations).

The algorithm still works, but we have to choose $m = 3$ columns out of $n = 5$ instead of $n = 2$ rows out of $m = 3$. There are $\binom{5}{3} = 10$ possibilities.

An implementation in Python looks like this:

First, draw $m$ numbers from $\{0, \dots, n-1\}$. Then solve for $\mathbf{x}$ via least squares. This ensures that we always find an $\mathbf{x}$ and we avoid checking for linear independence. The only numeric instability could come then from np.allclose(B @ x_B, b).

In the end, we obtain as before $\mathbf{c}^T\mathbf{x}^* = 6$ with $% $.

Both methods I presented here are of course not usable in normal problems. To solve real-world LPs, there exists for example the library CVXOPT.

Categories:

Updated: