Linear Programming and the Simplex Algorithm

In the last post in this series we saw some simple examples of linear programs, derived the concept of a dual linear program, and saw the duality theorem and the complementary slackness conditions which give a rough sketch of the stopping criterion for an algorithm. This time we’ll go ahead and write this algorithm for solving linear programs, and next time we’ll apply the algorithm to an industry-strength version of the nutrition problem we saw last time. The algorithm we’ll implement is called the simplex algorithm. It was the first algorithm for solving linear programs, invented in the 1940’s by George Dantzig, and it’s still the leading practical algorithm, and it was a key part of a Nobel Prize. It’s by far one of the most important algorithms ever devised.

As usual, we’ll post all of the code written in the making of this post on this blog’s Github page.

Slack variables and equality constraints

The simplex algorithm can solve any kind of linear program, but it only accepts a special form of the program as input. So first we have to do some manipulations. Recall that the primal form of a linear program was the following minimization problem.

\min \left \langle c, x \right \rangle \\ \textup{s.t. } Ax \geq b, x \geq 0

where the brackets mean “dot product.” And its dual is

\max \left \langle y, b \right \rangle \\ \textup{s.t. } A^Ty \leq c, y \geq 0

The linear program can actually have more complicated constraints than just the ones above. In general, one might want to have “greater than” and “less than” constraints in the same problem. It turns out that this isn’t any harder, and moreover the simplex algorithm only uses equality constraints, and with some finicky algebra we can turn any set of inequality or equality constraints into a set of equality constraints.

We’ll call our goal the “standard form,” which is as follows:

\max \left \langle c, x \right \rangle \\ \textup{s.t. } Ax = b, x \geq 0

It seems impossible to get the usual minimization/maximization problem into standard form until you realize there’s nothing stopping you from adding more variables to the problem. That is, say we’re given a constraint like:

\displaystyle x_7 + x_3 \leq 10,

we can add a new variable \xi, called a slack variable, so that we get an equality:

\displaystyle x_7 + x_3 + \xi = 10

And now we can just impose that \xi \geq 0. The idea is that \xi represents how much “slack” there is in the inequality, and you can always choose it to make the condition an equality. So if the equality holds and the variables are nonnegative, then the x_i will still satisfy their original inequality. For “greater than” constraints, we can do the same thing but subtract a nonnegative variable. Finally, if we have a minimization problem “\min z” we can convert it to \max -z.

So, to combine all of this together, if we have the following linear program with each kind of constraint,

Screen Shot 2014-10-05 at 12.06.19 AM

We can add new variables \xi_1, \xi_2, and write it as

Screen Shot 2014-10-05 at 12.06.41 AM

By defining the vector variable x = (x_1, x_2, x_3, \xi_1, \xi_2) and c = (-1,-1,-1,0,0) and A to have -1, 0, 1 as appropriately for the new variables, we see that the system is written in standard form.

This is the kind of tedious transformation we can automate with a program. Assuming there are n variables, the input consists of the vector c of length n, and three matrix-vector pairs (A, b) representing the three kinds of constraints. It’s a bit annoying to describe, but the essential idea is that we compute a rectangular “identity” matrix whose diagonal entries are \pm 1, and then join this with the original constraint matrix row-wise. The reader can see the full implementation in the Github repository for this post, though we won’t use this particular functionality in the algorithm that follows.

There are some other additional things we could do: for example there might be some variables that are completely unrestricted. What you do in this case is take an unrestricted variable z and replace it by the difference of two unrestricted variables z' - z''.  For simplicity we’ll ignore this, but it would be a fruitful exercise for the reader to augment the function to account for these.

What happened to the slackness conditions?

The “standard form” of our linear program raises an obvious question: how can the complementary slackness conditions make sense if everything is an equality? It turns out that one can redo all the work one did for linear programs of the form we gave last time (minimize w.r.t. greater-than constraints) for programs in the new “standard form” above. We even get the same complementary slackness conditions! If you want to, you can do this entire routine quite a bit faster if you invoke the power of Lagrangians. We won’t do that here, but the tool shows up as a way to work with primal-dual conversions in many other parts of mathematics, so it’s a good buzzword to keep in mind.

In our case, the only difference with the complementary slackness conditions is that one of the two is trivial: \left \langle y^*, Ax^* - b \right \rangle = 0. This is because if our candidate solution x^* is feasible, then it will have to satisfy Ax = b already. The other one, that \left \langle x^*, A^Ty^* - c \right \rangle = 0, is the only one we need to worry about.

Again, the complementary slackness conditions give us inspiration here. Recall that, informally, they say that when a variable is used at all, it is used as much as it can be to fulfill its constraint (the corresponding dual constraint is tight). So a solution will correspond to a choice of some variables which are either used or not, and a choice of nonzero variables will correspond to a solution. We even saw this happen in the last post when we observed that broccoli trumps oranges. If we can get a good handle on how to navigate the set of these solutions, then we’ll have a nifty algorithm.

Let’s make this official and lay out our assumptions.

Extreme points and basic solutions

Remember that the graphical way to solve a linear program is to look at the line (or hyperplane) given by \langle c, x \rangle = q and keep increasing q (or decreasing it, if you are minimizing) until the very last moment when this line touches the region of feasible solutions. Also recall that the “feasible region” is just the set of all solutions to Ax = b, that is the solutions that satisfy the constraints. We imagined this picture:

The constraints define a convex area of "feasible solutions." Image source: Wikipedia.

The constraints define a convex area of “feasible solutions.” Image source: Wikipedia.

With this geometric intuition it’s clear that there will always be an optimal solution on a vertex of the feasible region. These points are called extreme points of the feasible region. But because we will almost never work in the plane again (even introducing slack variables makes us relatively high dimensional!) we want an algebraic characterization of these extreme points.

If you have a little bit of practice with convex sets the correct definition is very natural. Recall that a set X is convex if for any two points x, y \in X every point on the line segment between x and y is also in X. An algebraic way to say this (thinking of these points now as vectors) is that every point \delta x + (1-\delta) y \in X when 0 \leq \delta \leq 1. Now an extreme point is just a point that isn’t on the inside of any such line, i.e. can’t be written this way for 0 < \delta < 1. For example,

A convex set with extremal points in red. Image credit Wikipedia.

A convex set with extremal points in red. Image credit Wikipedia.

Another way to say this is that if z is an extreme point then whenever z can be written as \delta x + (1-\delta) y for some 0 < \delta < 1, then actually x=y=z. Now since our constraints are all linear (and there are a finite number of them) they won’t define a convex set with weird curves like the one above. This means that there are a finite number of extreme points that just correspond to the intersections of some of the constraints. So there are at most 2^n possibilities.

Indeed we want a characterization of extreme points that’s specific to linear programs in standard form, “\max \langle c, x \rangle \textup{ s.t. } Ax=b, x \geq 0.” And here is one.

Definition: Let A be an m \times n matrix with n \geq m. A solution x to Ax=b is called basic if at most m of its entries are nonzero.

The reason we call it “basic” is because, under some mild assumptions we describe below, a basic solution corresponds to a vector space basis of \mathbb{R}^m. Which basis? The one given by the m columns of A used in the basic solution. We don’t need to talk about bases like this, though, so in the event of a headache just think of the basis as a set B \subset \{ 1, 2, \dots, n \} of size m corresponding to the nonzero entries of the basic solution.

Indeed, what we’re doing here is looking at the matrix A_B formed by taking the columns of A whose indices are in B, and the vector x_B in the same way, and looking at the equation A_Bx_B = b. If all the parts of x that we removed were zero then this will hold if and only if Ax=b. One might worry that A_B is not invertible, so we’ll go ahead and assume it is. In fact, we’ll assume that every set of m columns of A forms a basis and that the rows of A are also linearly independent. This isn’t without loss of generality because if some rows or columns are not linearly independent, we can remove the offending constraints and variables without changing the set of solutions (this is why it’s so nice to work with the standard form).

Moreover, we’ll assume that every basic solution has exactly m nonzero variables. A basic solution which doesn’t satisfy this assumption is called degenerate, and they’ll essentially be special corner cases in the simplex algorithm. Finally, we call a basic solution feasible if (in addition to satisfying Ax=b) it satisfies x \geq 0. Now that we’ve made all these assumptions it’s easy to see that choosing m nonzero variables uniquely determines a basic feasible solution. Again calling the sub-matrix A_B for a basis B, it’s just x_B = A_B^{-1}b. Now to finish our characterization, we just have to show that under the same assumptions basic feasible solutions are exactly the extremal points of the feasible region.

Proposition: A vector x is a basic feasible solution if and only if it’s an extreme point of the set \{ x : Ax = b, x \geq 0 \}.

Proof. For one direction, suppose you have a basic feasible solution x, and say we write it as x = \delta y + (1-\delta) z for some 0 < \delta < 1. We want to show that this implies y = z. Since all of these points are in the feasible region, all of their coordinates are nonnegative. So whenever a coordinate x_i = 0 it must be that both y_i = z_i = 0. Since x has exactly n-m zero entries, it must be that y, z both have at least n-m zero entries, and hence y,z are both basic. By our non-degeneracy assumption they both then have exactly m nonzero entries. Let B be the set of the nonzero indices of x. Because Ay = Az = b, we have A(y-z) = 0. Now y-z has all of its nonzero entries in B, and because the columns of A_B are linearly independent, the fact that A_B(y-z) = 0 implies y-z = 0.

In the other direction, suppose  that you have some extreme point x which is feasible but not basic. In other words, there are more than m nonzero entries of x, and we’ll call the indices J = \{ j_1, \dots, j_t \} where t > m. The columns of A_J are linearly dependent (since they’re t vectors in \mathbb{R}^m), and so let \sum_{i=1}^t z_{j_i} A_{j_i} be a nontrivial linear combination of the columns of A. Add zeros to make the z_{j_i} into a length n vector z, so that Az = 0. Now

A(x + \varepsilon z) = A(x - \varepsilon z) = Ax = b

And if we pick \varepsilon sufficiently small x \pm \varepsilon z will still be nonnegative, because the only entries we’re changing of x are the strictly positive ones. Then x = \delta (x + \varepsilon z) + (1 - \delta) \varepsilon z for \delta = 1/2, but this is very embarrassing for x who was supposed to be an extreme point. \square

Now that we know extreme points are the same as basic feasible solutions, we need to show that any linear program that has some solution has a basic feasible solution. This is clear geometrically: any time you have an optimum it has to either lie on a line or at a vertex, and if it lies on a line then you can slide it to a vertex without changing its value. Nevertheless, it is a useful exercise to go through the algebra.

Theorem. Whenever a linear program is feasible and bounded, it has a basic feasible solution.

Proof. Let x be an optimal solution to the LP. If x has at most m nonzero entries then it’s a basic solution and by the non-degeneracy assumption it must have exactly m nonzero entries. In this case there’s nothing to do, so suppose that x has r > m nonzero entries. It can’t be a basic feasible solution, and hence is not an extreme point of the set of feasible solutions (as proved by the last theorem). So write it as x = \delta y + (1-\delta) z for some feasible y \neq z and 0 < \delta < 1.

The only thing we know about x is it’s optimal. Let c be the cost vector, and the optimality says that \langle c,x \rangle \geq \langle c,y \rangle, and \langle c,x \rangle \geq \langle c,z \rangle. We claim that in fact these are equal, that y, z are both optimal as well. Indeed, say y were not optimal, then

\displaystyle \langle c, y \rangle < \langle c,x \rangle = \delta \langle c,y \rangle + (1-\delta) \langle c,z \rangle

Which can be rearranged to show that \langle c,y \rangle < \langle c, z \rangle. Unfortunately for x, this implies that it was not optimal all along:

\displaystyle \langle c,x \rangle < \delta \langle c, z \rangle + (1-\delta) \langle c,z \rangle = \langle c,z \rangle

An identical argument works to show z is optimal, too. Now we claim we can use y,z to get a new solution that has fewer than r nonzero entries. Once we show this we’re done: inductively repeat the argument with the smaller solution until we get down to exactly m nonzero variables. As before we know that y,z must have at least as many zeros as x. If they have more zeros we’re done. And if they have exactly as many zeros we can do the following trick. Write w = \gamma y + (1- \gamma)z for a \gamma \in \mathbb{R} we’ll choose later. Note that no matter the \gamma, w is optimal. Rewriting w = z + \gamma (y-z), we just have to pick a \gamma that ensures one of the nonzero coefficients of z is zeroed out while maintaining nonnegativity. Indeed, we can just look at the index i which minimizes z_i / (y-z)_i and use \delta = - z_i / (y-z)_i. \square.

So we have an immediate (and inefficient) combinatorial algorithm: enumerate all subsets of size m, compute the corresponding basic feasible solution x_B = A_B^{-1}b, and see which gives the biggest objective value. The problem is that, even if we knew the value of m, this would take time n^m, and it’s not uncommon for m to be in the tens or hundreds (and if we don’t know m the trivial search is exponential).

So we have to be smarter, and this is where the simplex tableau comes in.

The simplex tableau

Now say you have any basis B and any feasible solution x. For now x might not be a basic solution, and even if it is, its basis of nonzero entries might not be the same as B. We can decompose the equation Ax = b into the basis part and the non basis part:

A_Bx_B + A_{B'} x_{B'} = b

and solving the equation for x_B gives

x_B = A^{-1}_B(b - A_{B'} x_{B'})

It may look like we’re making a wicked abuse of notation here, but both A_Bx_B and A_{B'}x_{B'} are vectors of length m so the dimensions actually do work out. Now our feasible solution x has to satisfy Ax = b, and the entries of x are all nonnegative, so it must be that x_B \geq 0 and x_{B'} \geq 0, and by the equality above A^{-1}_B (b - A_{B'}x_{B'}) \geq 0 as well. Now let’s write the maximization objective \langle c, x \rangle by expanding it first in terms of the x_B, x_{B'}, and then expanding x_B.

\displaystyle \begin{aligned} \langle c, x \rangle & = \langle c_B, x_B \rangle + \langle c_{B'}, x_{B'} \rangle \\  & = \langle c_B, A^{-1}_B(b - A_{B'}x_{B'}) \rangle + \langle c_{B'}, x_{B'} \rangle \\  & = \langle c_B, A^{-1}_Bb \rangle + \langle c_{B'} - (A^{-1}_B A_{B'})^T c_B, x_{B'} \rangle \end{aligned}

If we want to maximize the objective, we can just maximize this last line. There are two cases. In the first, the vector c_{B'} - (A^{-1}_B A_{B'})^T c_B \leq 0 and A_B^{-1}b \geq 0. In the above equation, this tells us that making any component of x_{B'} bigger will decrease the overall objective. In other words, \langle c, x \rangle \leq \langle c_B, A_B^{-1}b \rangle. Picking x = A_B^{-1}b (with zeros in the non basis part) meets this bound and hence must be optimal. In other words, no matter what basis B we’ve chosen (i.e., no matter the candidate basic feasible solution), if the two conditions hold then we’re done.

Now the crux of the algorithm is the second case: if the conditions aren’t met, we can pick a positive index of c_{B'} - (A_B^{-1}A_{B'})^Tc_B and increase the corresponding value of x_{B'} to increase the objective value. As we do this, other variables in the solution will change as well (by decreasing), and we have to stop when one of them hits zero. In doing so, this changes the basis by removing one index and adding another. In reality, we’ll figure out how much to increase ahead of time, and the change will correspond to a single elementary row-operation in a matrix.

Indeed, the matrix we’ll use to represent all of this data is called a tableau in the literature. The columns of the tableau will correspond to variables, and the rows to constraints. The last row of the tableau will maintain a candidate solution y to the dual problem. Here’s a rough picture to keep the different parts clear while we go through the details.

tableau

But to make it work we do a slick trick, which is to “left-multiply everything” by A_B^{-1}. In particular, if we have an LP given by c, A, b, then for any basis it’s equivalent to the LP given by c, A_B^{-1}A, A_{B}^{-1} b (just multiply your solution to the new program by A_B to get a solution to the old one). And so the actual tableau will be of this form.

tableau-symbols

When we say it’s in this form, it’s really only true up to rearranging columns. This is because the chosen basis will always be represented by an identity matrix (as it is to start with), so to find the basis you can find the embedded identity sub-matrix. In fact, the beginning of the simplex algorithm will have the initial basis sitting in the last few columns of the tableau.

Let’s look a little bit closer at the last row. The first portion is zero because A_B^{-1}A_B is the identity. But furthermore with this A_B^{-1} trick the dual LP involves A_B^{-1} everywhere there’s a variable. In particular, joining all but the last column of the last row of the tableau, we have the vector c - A_B^T(A_B^{-1})^T c, and setting y = A_B^{-1}c_B we get a candidate solution for the dual. What makes the trick even slicker is that A_B^{-1}b is already the candidate solution x_B, since (A_B^{-1}A)_B^{-1} is the identity. So we’re implicitly keeping track of two solutions here, one for the primal LP, given by the last column of the tableau, and one for the dual, contained in the last row of the tableau.

I told you the last row was the dual solution, so why all the other crap there? This is the final slick in the trick: the last row further encodes the complementary slackness conditions. Now that we recognize the dual candidate sitting there, the complementary slackness conditions simply ask for the last row to be non-positive (this is just another way of saying what we said at the beginning of this section!). You should check this, but it gives us a stopping criterion: if the last row is non-positive then stop and output the last column.

The simplex algorithm

Now (finally!) we can describe and implement the simplex algorithm in its full glory. Recall that our informal setup has been:

  1. Find an initial basic feasible solution, and set up the corresponding tableau.
  2. Find a positive index of the last row, and increase the corresponding variable (adding it to the basis) just enough to make another variable from the basis zero (removing it from the basis).
  3. Repeat step 2 until the last row is nonpositive.
  4. Output the last column.

This is almost correct, except for some details about how increasing the corresponding variables works. What we’ll really do is represent the basis variables as pivots (ones in the tableau) and then the first 1 in each row will be the variable whose value is given by the entry in the last column of that row. So, for example, the last entry in the first row may be the optimal value for x_5, if the fifth column is the first entry in row 1 to have a 1.

As we describe the algorithm, we’ll illustrate it running on a simple example. In doing this we’ll see what all the different parts of the tableau correspond to from the previous section in each step of the algorithm.

example

Spoiler alert: the optimum is x_1 = 2, x_2 = 1 and the value of the max is 8.

So let’s be more programmatically formal about this. The main routine is essentially pseudocode, and the difficulty is in implementing the helper functions

def simplex(c, A, b):
   tableau = initialTableau(c, A, b)

   while canImprove(tableau):
      pivot = findPivotIndex(tableau)
      pivotAbout(tableau, pivot)

   return primalSolution(tableau), objectiveValue(tableau)

Let’s start with the initial tableau. We’ll assume the user’s inputs already include the slack variables. In particular, our example data before adding slack is

c = [3, 2]
A = [[1, 2], [1, -1]]
b = [4, 1]

And after adding slack:

c = [3, 2, 0, 0]
A = [[1,  2,  1,  0],
     [1, -1,  0,  1]]
b = [4, 1]

Now to set up the initial tableau we need an initial feasible solution in mind. The reader is recommended to work this part out with a pencil, since it’s much easier to write down than it is to explain. Since we introduced slack variables, our initial feasible solution (basis) B can just be (0,0,1,1). And so x_B is just the slack variables, c_B is the zero vector, and A_B is the 2×2 identity matrix. Now A_B^{-1}A_{B'} = A_{B'}, which is just the original two columns of A we started with, and A_B^{-1}b = b. For the last row, c_B is zero so the part under A_B^{-1}A_B is the zero vector. The part under A_B^{-1}A_{B'} is just c_{B'} = (3,2).

Rather than move columns around every time the basis B changes, we’ll keep the tableau columns in order of (x_1, \dots, x_n, \xi_1, \dots, \xi_m). In other words, for our example the initial tableau should look like this.

[[ 1,  2,  1,  0,  4],
 [ 1, -1,  0,  1,  1],
 [ 3,  2,  0,  0,  0]]

So implementing initialTableau is just a matter of putting the data in the right place.

def initialTableau(c, A, b):
   tableau = [row[:] + [x] for row, x in zip(A, b)]
   tableau.append(c[:] + [0])
   return tableau

As an aside: in the event that we don’t start with the trivial basic feasible solution of “trivially use the slack variables,” we’d have to do a lot more work in this function. Next, the primalSolution() and objectiveValue() functions are simple, because they just extract the encoded information out from the tableau (some helper functions are omitted for brevity).

def primalSolution(tableau):
   # the pivot columns denote which variables are used
   columns = transpose(tableau)
   indices = [j for j, col in enumerate(columns[:-1]) if isPivotCol(col)]
   return list(zip(indices, columns[-1]))

def objectiveValue(tableau):
   return -(tableau[-1][-1])

Similarly, the canImprove() function just checks if there’s a nonnegative entry in the last row

def canImprove(tableau):
   lastRow = tableau[-1]
   return any(x > 0 for x in lastRow[:-1])

Let’s run the first loop of our simplex algorithm. The first step is checking to see if anything can be improved (in our example it can). Then we have to find a pivot entry in the tableau. This part includes some edge-case checking, but if the edge cases aren’t a problem then the strategy is simple: find a positive entry corresponding to some entry j of B', and then pick an appropriate entry in that column to use as the pivot. Pivoting increases the value of x_j (from zero) to whatever is the largest we can make it without making some other variables become negative. As we’ve said before, we’ll stop increasing x_j when some other variable hits zero, and we can compute which will be the first to do so by looking at the current values of x_B = A_B^{-1}b (in the last column of the tableau), and seeing how pivoting will affect them. If you stare at it for long enough, it becomes clear that the first variable to hit zero will be the entry x_i of the basis for which x_i / A_{i,j} is minimal (and A_{i,j} has to be positve). This is because, in order to maintain the linear equalities, every entry of x_B will be decreased by that value during a pivot, and we can’t let any of the variables become negative.

All of this results in the following function, where we have left out the degeneracy/unboundedness checks.

def findPivotIndex(tableau):
   # pick first nonzero index of the last row
   column = [i for i,x in enumerate(tableau[-1][:-1]) if x > 0][0]
   quotients = [(i, r[-1] / r[column]) for i,r in enumerate(tableau[:-1]) if r[column] > 0]

   # pick row index minimizing the quotient
   row = min(quotients, key=lambda x: x[1])[0]
   return row, column

For our example, the minimizer is the (1,0) entry (second row, first column). Pivoting is just doing the usual elementary row operations (we covered this in a primer a while back on row-reduction). The pivot function we use here is no different, and in particular mutates the list in place.

def pivotAbout(tableau, pivot):
   i,j = pivot

   pivotDenom = tableau[i][j]
   tableau[i] = [x / pivotDenom for x in tableau[i]]

   for k,row in enumerate(tableau):
      if k != i:
         pivotRowMultiple = [y * tableau[k][j] for y in tableau[i]]
         tableau[k] = [x - y for x,y in zip(tableau[k], pivotRowMultiple)]

And in our example pivoting around the chosen entry gives the new tableau.

[[ 0.,  3.,  1., -1.,  3.],
 [ 1., -1.,  0.,  1.,  1.],
 [ 0.,  5.,  0., -3., -3.]]

In particular, B is now (1,0,1,0), since our pivot removed the second slack variable \xi_2 from the basis. Currently our solution has x_1 = 1, \xi_1 = 3. Notice how the identity submatrix is still sitting in there, the columns are just swapped around.

There’s still a positive entry in the bottom row, so let’s continue. The next pivot is (0,1), and pivoting around that entry gives the following tableau:

[[ 0.        ,  1.        ,  0.33333333, -0.33333333,  1.        ],
 [ 1.        ,  0.        ,  0.33333333,  0.66666667,  2.        ],
 [ 0.        ,  0.        , -1.66666667, -1.33333333, -8.        ]]

And because all of the entries in the bottom row are negative, we’re done. We read off the solution as we described, so that the first variable is 2 and the second is 1, and the objective value is the opposite of the bottom right entry, 8.

To see all of the source code, including the edge-case-checking we left out of this post, see the Github repository for this post.

Obivous questions and sad answers

An obvious question is: what is the runtime of the simplex algorithm? Is it polynomial in the size of the tableau? Is it even guaranteed to stop at some point? The surprising truth is that nobody knows the answer to all of these questions! Originally (in the 1940’s) the simplex algorithm actually had an exponential runtime in the worst case, though this was not known until 1972. And indeed, to this day while some variations are known to terminate, no variation is known to have polynomial runtime in the worst case. Some of the choices we made in our implementation (for example, picking the first column with a positive entry in the bottom row) have the potential to cycle, i.e., variables leave and enter the basis without changing the objective at all. Doing something like picking a random positive column, or picking the column which will increase the objective value by the largest amount are alternatives. Unfortunately, every single pivot-picking rule is known to give rise to exponential-time simplex algorithms in the worst case (in fact, this was discovered as recently as 2011!). So it remains open whether there is a variant of the simplex method that runs in guaranteed polynomial time.

But then, in a stunning turn of events, Leonid Khachiyan proved in the 70’s that in fact linear programs can always be solved in polynomial time, via a completely different algorithm called the ellipsoid method. Following that was a method called the interior point method, which is significantly more efficient. Both of these algorithms generalize to problems that are harder than linear programming as well, so we will probably cover them in the distant future of this blog.

Despite the celebratory nature of these two results, people still use the simplex algorithm for industrial applications of linear programming. The reason is that it’s much faster in practice, and much simpler to implement and experiment with.

The next obvious question has to do with the poignant observation that whole numbers are great. That is, you often want the solution to your problem to involve integers, and not real numbers. But adding the constraint that the variables in a linear program need to be integer valued (even just 0-1 valued!) is NP-complete. This problem is called integer linear programming, or just integer programming (IP). So we can’t hope to solve IP, and rightly so: the reader can verify easily that boolean satisfiability instances can be written as linear programs where each clause corresponds to a constraint.

This brings up a very interesting theoretical issue: if we take an integer program and just remove the integrality constraints, and solve the resulting linear program, how far away are the two solutions? If they’re close, then we can hope to give a good approximation to the integer program by solving the linear program and somehow turning the resulting solution back into an integer solution. In fact this is a very popular technique called LP-rounding. We’ll also likely cover that on this blog at some point.

Oh there’s so much to do and so little time! Until next time.

About these ads

Fixing Bugs in “Computing Homology”

A few awesome readers have posted comments in Computing Homology to the effect of, “Your code is not quite correct!” And they’re right! Despite the almost year since that post’s publication, I haven’t bothered to test it for more complicated simplicial complexes, or even the basic edge cases! When I posted it the mathematics just felt so solid to me that it had to be right (the irony is rich, I know).

As such I’m apologizing for my lack of rigor and explaining what went wrong, the fix, and giving some test cases. As of the publishing of this post, the Github repository for Computing Homology has been updated with the correct code, and some more examples.

The main subroutine was the simultaneousReduce function which I’ll post in its incorrectness below

def simultaneousReduce(A, B):
   if A.shape[1] != B.shape[0]:
      raise Exception("Matrices have the wrong shape.")

   numRows, numCols = A.shape # col reduce A

   i,j = 0,0
   while True:
      if i >= numRows or j >= numCols:
         break

      if A[i][j] == 0:
         nonzeroCol = j
         while nonzeroCol < numCols and A[i,nonzeroCol] == 0:
            nonzeroCol += 1

         if nonzeroCol == numCols:
            j += 1
            continue

         colSwap(A, j, nonzeroCol)
         rowSwap(B, j, nonzeroCol)

      pivot = A[i,j]
      scaleCol(A, j, 1.0 / pivot)
      scaleRow(B, j, 1.0 / pivot)

      for otherCol in range(0, numCols):
         if otherCol == j:
            continue
         if A[i, otherCol] != 0:
            scaleAmt = -A[i, otherCol]
            colCombine(A, otherCol, j, scaleAmt)
            rowCombine(B, j, otherCol, -scaleAmt)

      i += 1; j+= 1

   return A,B

It’s a beast of a function, and the persnickety detail was just as beastly: this snippet should have an i += 1 instead of a j.

if nonzeroCol == numCols:
   j += 1
   continue

This is simply what happens when we’re looking for a nonzero entry in a row to use as a pivot for the corresponding column, but we can’t find one and have to move to the next row. A stupid error on my part that would be easily caught by proper test cases.

The next mistake is a mathematical misunderstanding. In short, the simultaneous column/row reduction process is not enough to get the \partial_{k+1} matrix into the right form! Let’s see this with a nice example, a triangulation of the Mobius band. There are a number of triangulations we could use, many of which are seen in these slides. The one we’ll use is the following.

mobius-triangulation

It’s first and second boundary maps are as follows (in code, because latex takes too much time to type out)

mobiusD1 = numpy.array([
   [-1,-1,-1,-1, 0, 0, 0, 0, 0, 0],
   [ 1, 0, 0, 0,-1,-1,-1, 0, 0, 0],
   [ 0, 1, 0, 0, 1, 0, 0,-1,-1, 0],
   [ 0, 0, 0, 1, 0, 0, 1, 0, 1, 1],
])

mobiusD2 = numpy.array([
   [ 1, 0, 0, 0, 1],
   [ 0, 0, 0, 1, 0],
   [-1, 0, 0, 0, 0],
   [ 0, 0, 0,-1,-1],
   [ 0, 1, 0, 0, 0],
   [ 1,-1, 0, 0, 0],
   [ 0, 0, 0, 0, 1],
   [ 0, 1, 1, 0, 0],
   [ 0, 0,-1, 1, 0],
   [ 0, 0, 1, 0, 0],
])

And if we were to run the above code on it we’d get a first Betti number of zero (which is incorrect, it’s first homology group has rank 1). Here are the reduced matrices.

>>> A1, B1 = simultaneousReduce(mobiusD1, mobiusD2)
>>> A1
array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]])
>>> B1
array([[ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  1,  0,  0,  0],
       [ 1, -1,  0,  0,  0],
       [ 0,  0,  0,  0,  1],
       [ 0,  1,  1,  0,  0],
       [ 0,  0, -1,  1,  0],
       [ 0,  0,  1,  0,  0]])

The first reduced matrix looks fine; there’s nothing we can do to improve it. But the second one is not quite fully reduced! Notice that rows 5, 8 and 10 are not linearly independent. So we need to further row-reduce the nonzero part of this matrix before we can read off the true rank in the way we described last time. This isn’t so hard (we just need to reuse the old row-reduce function we’ve been using), but why is this allowed? It’s just because the corresponding column operations for those row operations are operating on columns of all zeros! So we need not worry about screwing up the work we did in column reducing the first matrix, as long as we only work with the nonzero rows of the second.

Of course, nothing is stopping us from ignoring the “corresponding” column operations, since we know we’re already done there. So we just have to finish row reducing this matrix.

This changes our bettiNumber function by adding a single call to a row-reduce function which we name so as to be clear what’s happening. The resulting function is

def bettiNumber(d_k, d_kplus1):
   A, B = numpy.copy(d_k), numpy.copy(d_kplus1)
   simultaneousReduce(A, B)
   finishRowReducing(B)

   dimKChains = A.shape[1]
   kernelDim = dimKChains - numPivotCols(A)
   imageDim = numPivotRows(B)

   return kernelDim - imageDim

And running this on our Mobius band example gives:

>>> bettiNumber(mobiusD1, mobiusD2))
1

As desired. Just to make sure things are going swimmingly under the hood, we can check to see how finishRowReducing does after calling simultaneousReduce

>>> simultaneousReduce(mobiusD1, mobiusD2)
(array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]]), array([[ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  1,  0,  0,  0],
       [ 1, -1,  0,  0,  0],
       [ 0,  0,  0,  0,  1],
       [ 0,  1,  1,  0,  0],
       [ 0,  0, -1,  1,  0],
       [ 0,  0,  1,  0,  0]]))
>>> finishRowReducing(mobiusD2)
array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0]])

Indeed, finishRowReducing finishes row reducing the second boundary matrix. Note that it doesn’t preserve how the rows of zeros lined up with the pivot columns of the reduced version of \partial_1 as it did in the previous post, but since in the end we’re only counting pivots it doesn’t matter how we switch rows. The “zeros lining up” part is just for a conceptual understanding of how the image lines up with the kernel for a valid simplicial complex.

In fixing this issue we’ve also fixed an issue another commenter mentioned, that you couldn’t blindly plug in the zero matrix for \partial_0 and get zeroth homology (which is the same thing as connected components). After our fix you can.

Of course there still might be bugs, but I have so many drafts lined up on this blog (and research papers to write, experiments to run, theorems to prove), that I’m going to put off writing a full test suite. I’ll just have to update this post with new bug fixes as they come. There’s just so much math and so little time :) But extra kudos to my amazing readers who were diligent enough to run examples and spot my error. I’m truly blessed to have you on my side.

Also note that this isn’t the most efficient way to represent the simplicial complex data, or the most efficient row reduction algorithm. If you’re going to run the code on big inputs, I suggest you take advantage of sparse matrix algorithms for doing this sort of stuff. You can represent the simplices as entries in a dictionary and do all sorts of clever optimizations to make the algorithm effectively linear time in the number of simplices.

Until next time!

Computing Homology

Update: the mistakes made in the code posted here are fixed and explained in a subsequent post (one minor code bug was fixed here, and a less minor conceptual bug is fixed in the linked post).

In our last post in this series on topology, we defined the homology group. Specifically, we built up a topological space as a simplicial complex (a mess of triangles glued together), we defined an algebraic way to represent collections of simplices called chains as vectors in a vector space, we defined the boundary homomorphism \partial_k as a linear map on chains, and finally defined the homology groups as the quotient vector spaces

\displaystyle H_k(X) = \frac{\textup{ker} \partial_k}{\textup{im} \partial_{k+1}}.

The number of holes in X was just the dimension of this quotient space.

In this post we will be quite a bit more explicit. Because the chain groups are vector spaces and the boundary mappings are linear maps, they can be represented as matrices whose dimensions depend on our simplicial complex structure. Better yet, if we have explicit representations of our chains by way of a basis, then we can use row-reduction techniques to write the matrix in a standard form.

Of course the problem arises when we want to work with two matrices simultaneously (to compute the kernel-mod-image quotient above). This is not computationally any more difficult, but it requires some theoretical fiddling. We will need to dip a bit deeper into our linear algebra toolboxes to see how it works, so the rusty reader should brush up on their linear algebra before continuing (or at least take some time to sort things out if or when confusion strikes).

Without further ado, let’s do an extended example and work our ways toward a general algorithm. As usual, all of the code used for this post is available on this blog’s Github page.

Two Big Matrices

Recall our example simplicial complex from last time.

circle-wedge-sphere

We will compute H_1 of this simplex (which we saw last time was \mathbb{Q}) in a more algorithmic way than we did last time.

Once again, we label the vertices 0-4 so that the extra “arm” has vertex 4 in the middle, and its two endpoints are 0 and 2. This gave us orientations on all of the simplices, and the following chain groups. Since the vertex labels (and ordering) are part of the data of a simplicial complex, we have made no choices in writing these down.

\displaystyle C_0(X) = \textup{span} \left \{ 0,1,2,3,4 \right \}

\displaystyle C_1(X) = \textup{span} \left \{ [0,1], [0,2], [0,3], [0,4], [1,2], [1,3],[2,3],[2,4] \right \}

\displaystyle C_2(X) = \textup{span} \left \{ [0,1,2], [0,1,3], [0,2,3], [1,2,3] \right \}

Now given our known definitions of \partial_k as an alternating sum from last time, we can give a complete specification of the boundary map as a matrix. For \partial_1, this would be

\displaystyle \partial_1 = \bordermatrix{  & [0,1] & [0,2] & [0,3] & [0,4] & [1,2] & [1,3] & [2,3] & [2,4] \cr  0 & -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0\cr  1 & 1 & 0 & 0 & 0 & -1 & -1 & 0 & 0\cr  2 & 0 & 1 & 0 & 0 & 1 & 0 & -1 & -1 \cr  3 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 \cr  4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 },

where the row labels are the basis for C_0(X) and the column labels are the basis for C_1(X). Similarly, \partial_2 is

\displaystyle \partial_2 = \bordermatrix{  & [0,1,2] & [0,1,3] & [0,2,3] & [1,2,3] \cr  [0,1] & 1 & 1 & 0 & 0\cr  [0,2] & -1 & 0 & 1 & 0\cr  [0,3] & 0 & -1 & -1 & 0\cr  [0,4] & 0 & 0 & 0 & 0\cr  [1,2] & 1 & 0 & 0 & 1\cr  [1,3] & 0 & 1 & 0 & -1\cr  [2,3] & 0 & 0 & 1 & 1\cr  [2,4] & 0 & 0 & 0 & 0}

The reader is encouraged to check that these matrices are written correctly by referring to the formula for \partial as given last time.

Remember the crucial property of \partial, that \partial^2 = \partial_k \partial_{k+1} = 0. Indeed, the composition of the two boundary maps just corresponds to the matrix product of the two matrices, and one can verify by hand that the above two matrices multiply to the zero matrix.

We know from basic linear algebra how to compute the kernel of a linear map expressed as a matrix: column reduce and inspect the columns of zeros. Since the process of row reducing is really a change of basis, we can encapsulate the reduction inside a single invertible matrix A, which, when left-multiplied by \partial, gives us the reduced form of the latter. So write the reduced form of \partial_1 as \partial_1 A.

However, now we’re using two different sets of bases for the shared vector space involved in \partial_1 and \partial_2. In general, it will no longer be the case that \partial_kA\partial_{k+1} = 0. The way to alleviate this is to perform the “corresponding” change of basis in \partial_{k+1}. To make this idea more transparent, we return to the basics.

Changing Two Matrices Simultaneously

Recall that a matrix M represents a linear map between two vector spaces f : V \to W. The actual entries of M depend crucially on the choice of a basis for the domain and codomain. Indeed, if v_i form a basis for V and w_j for W, then the k-th column of the matrix representation M is defined to be the coefficients of the representation of f(v_k) in terms of the w_j. We hope to have nailed this concept down firmly in our first linear algebra primer.

Recall further that row operations correspond to changing a basis for the codomain, and column operations correspond to changing a basis for the domain. For example, the idea of swapping columns i,j in M gives a new matrix which is the representation of f with respect to the (ordered) basis for V which swaps the order of v_i , v_j. Similar things happen for all column operations (they all correspond to manipulations of the basis for V), while analogously row operations implicitly transform the basis for the codomain. Note, though, that the connection between row operations and transformations of the basis for the codomain are slightly more complicated than they are for the column operations. We will explicitly see how it works later in the post.

And so if we’re working with two maps A: U \to V and B: V \to W, and we change a basis for V in B via column reductions, then in order to be consistent, we need to change the basis for V in A via “complementary” row reductions. That is, if we call the change of basis matrix Q, then we’re implicitly sticking Q in between the composition BA to get (BQ)A. This is not the same map as BA, but we can make it the same map by adding a Q^{-1} in the right place:

\displaystyle BA = B(QQ^{-1})A = (BQ)(Q^{-1}A)

Indeed, whenever Q is a change of basis matrix so is Q^{-1} (trivially), and moreover the operations that Q performs on the columns of B are precisely the operations that Q^{-1} performs on the rows of A (this is because elementary row operations take different forms when multiplied on the left or right).

Coming back to our boundary operators, we want a canonical way to view the image of \partial_{k+1} as sitting inside the kernel of \partial_k. If we go ahead and use column reductions to transform \partial_k into a form where the kernel is easy to read off (as the columns consisting entirely of zeroes), then the corresponding row operations, when performed on \partial_{k+1} will tell us exactly the image of \partial_{k+1} inside the kernel of \partial_k.

This last point is true precisely because \textup{im} \partial_{k+1} \subset \textup{ker} \partial_k. This fact guarantees that the irrelevant rows of the reduced version of \partial_{k+1} are all zero.

Let’s go ahead and see this in action on our two big matrices above. For \partial_1, the column reduction matrix is

\displaystyle A =  \begin{pmatrix}  0 & 1 & 0 & 0 & 1 & 1 & 0 & 0\\  0 & 0 & 1 & 0 & -1 & 0 & 1 & 1\\  0 & 0 & 0 & 1 & 0 & -1 & -1 & 0\\  -1 & -1 & -1 & -1 & 0 & 0 & 0 & -1\\  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1  \end{pmatrix}

And the product \partial_1 A is

\displaystyle \partial_1 A =  \begin{pmatrix}  1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\  0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\  0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\  -1 & -1 & -1 & -1 & 0 & 0 & 0 & 0  \end{pmatrix}

Now the inverse of A, which is the corresponding basis change for \partial_2, is

\displaystyle A^{-1} =  \begin{pmatrix}  -1 & -1 & -1 & -1 & -0 & -0 & -0 & -0\\  1 & 0 & 0 & 0 & -1 & -1 & 0 & 0\\  0 & 1 & 0 & 0 & 1 & 0 & -1 & -1\\  0 & 0 & 1 & 0 & 0 & 1 & 1 & 0\\  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1  \end{pmatrix}

and the corresponding reduced form of \partial_2 is

\displaystyle A^{-1} \partial_2 =  \begin{pmatrix}  0 & 0 & 0 & 0\\  0 & 0 & 0 & 0\\  0 & 0 & 0 & 0\\  0 & 0 & 0 & 0\\  1 & 0 & 0 & 1\\  0 & 1 & 0 & -1\\  0 & 0 & 1 & 1\\  0 & 0 & 0 & 0  \end{pmatrix}

As a side note, we got these matrices by slightly modifying the code from our original post on row reduction to output the change of basis matrix in addition to performing row reduction. It turns out one can implement column reduction as row reduction of the transpose, and the change of basis matrix you get from this process will be the transpose of the change of basis matrix you want (by (AB)^\textup{T} = (B^\textup{T}A^\textup{T})). Though the code is particularly ad-hoc, we include it with the rest of the code used in this post on this blog’s Github page.

Now let’s inspect the two matrices \partial_1 A and A^{-1} \partial_2 more closely. The former has four “pivots” left over, and this corresponds to the rank of the matrix being 4. Moreover, the four basis vectors representing the columns with nonzero pivots, which we’ll call v_1, v_2, v_3, v_4 (we don’t care what their values are), span a complementary subspace to the kernel of \partial_1. Hence, the remaining four vectors (which we’ll call v_5, v_6, v_7, v_8) span the kernel. In particular, this says that the kernel has dimension 4.

On the other hand, we performed the same transformation of the basis of C_1(X) for \partial_2. Looking at the matrix that resulted, we see that the first four rows and the last row (representing v_1, v_2, v_3, v_4, v_8) are entirely zeros and so the image of \partial_2 intersects their span trivially. and the remaining three rows (representing v_5, v_6, v_7) have nonzero pivots. This tells us exactly that the image of \partial_2 is spanned by v_5, v_6, v_7.

And now, the coup de grâce, the quotient to get homology is simply

\displaystyle \frac{ \textup{span} \left \{ v_5, v_6, v_7, v_8 \right \}}{ \textup{span} \left \{ v_5, v_6, v_7 \right \}} = \textup{span} \left \{ v_8 \right \}

And the dimension of the homology group is 1, as desired.

The General Algorithm

It is no coincidence that things worked out at nicely as they did. The process we took of simultaneously rewriting two matrices with respect to a common basis is the bulk of the algorithm to compute homology. Since we’re really only interested in the dimensions of the homology groups, we just need to count pivots. If the number of pivots arising in \partial_k is y and the number of pivots arising in \partial_{k+1} is z, and the dimension of C_k(X) is n, then the dimension is exactly

(n-y) - z = \textup{dim}(\textup{ker} \partial_k) - \textup{dim}(\textup{im}\partial_{k+1})

And it is no coincidence that the pivots lined up so nicely to allow us to count dimensions this way. It is a minor exercise to prove it formally, but the fact that the composition \partial_k \partial_{k+1} = 0 implies that the reduced version of \partial_{k+1} will have an almost reduced row-echelon form (the only difference being the rows of zeros interspersed above, below, and possibly between pivot rows).

As the reader may have guessed at this point, we don’t actually need to compute A and A^{-1}. Instead of this, we can perform the column/row reductions simultaneously on the two matrices. The above analysis helped us prove the algorithm works, and with that guarantee we can throw out the analytical baggage and just compute the damn thing.

Indeed, assuming the input is already processed as two matrices representing the boundary operators with respect to the standard bases of the chain groups, computing homology is only slightly more difficult than row reducing in the first place. Putting our homology where our mouth is, we’ve implemented the algorithm in Python. As usual, the entire code used in this post is available on this blog’s Github page.

The first step is writing auxiliary functions to do elementary row and column operations on matrices. For this post, we will do everything in numpy (which makes the syntax shorter than standard Python syntax, but dependent on the numpy library).

import numpy

def rowSwap(A, i, j):
   temp = numpy.copy(A[i, :])
   A[i, :] = A[j, :]
   A[j, :] = temp

def colSwap(A, i, j):
   temp = numpy.copy(A[:, i])
   A[:, i] = A[:, j]
   A[:, j] = temp

def scaleCol(A, i, c):
   A[:, i] *= c*numpy.ones(A.shape[0])

def scaleRow(A, i, c):
   A[i, :] *= c*numpy.ones(A.shape[1])

def colCombine(A, addTo, scaleCol, scaleAmt):
   A[:, addTo] += scaleAmt * A[:, scaleCol]

def rowCombine(A, addTo, scaleRow, scaleAmt):
   A[addTo, :] += scaleAmt * A[scaleRow, :]

From here, the main meat of the algorithm is doing column reduction on one matrix, and applying the corresponding row operations on the other.

def simultaneousReduce(A, B):
   if A.shape[1] != B.shape[0]:
      raise Exception("Matrices have the wrong shape.")

   numRows, numCols = A.shape # col reduce A

   i,j = 0,0
   while True:
      if i >= numRows or j >= numCols:
         break

      if A[i][j] == 0:
         nonzeroCol = j
         while nonzeroCol < numCols and A[i,nonzeroCol] == 0:
            nonzeroCol += 1

         if nonzeroCol == numCols:
            i += 1
            continue

         colSwap(A, j, nonzeroCol)
         rowSwap(B, j, nonzeroCol)

      pivot = A[i,j]
      scaleCol(A, j, 1.0 / pivot)
      scaleRow(B, j, 1.0 / pivot)

      for otherCol in range(0, numCols):
         if otherCol == j:
            continue
         if A[i, otherCol] != 0:
            scaleAmt = -A[i, otherCol]
            colCombine(A, otherCol, j, scaleAmt)
            rowCombine(B, j, otherCol, -scaleAmt)

      i += 1; j+= 1

   return A,B

This more or less parallels the standard algorithm for row-reduction (with the caveat that all the indices are swapped to do column-reduction). The only somewhat confusing line is the call to rowCombine, which explicitly realizes the corresponding row operation as the inverse of the performed column operation. Note that for row operations, the correspondence between operations on the basis and operations on the rows is not as direct as it is for columns. What’s given above is the true correspondence. Writing down lots of examples will reveal why, and we leave that as an exercise to the reader.

Then the actual algorithm to compute homology is just a matter of counting pivots. Here are two pivot-counting functions in a typical numpy fashion:

def numPivotCols(A):
   z = numpy.zeros(A.shape[0])
   return [numpy.all(A[:, j] == z) for j in range(A.shape[1])].count(False)

def numPivotRows(A):
   z = numpy.zeros(A.shape[1])
   return [numpy.all(A[i, :] == z) for i in range(A.shape[0])].count(False)

And the final function is just:

def bettiNumber(d_k, d_kplus1):
   A, B = numpy.copy(d_k), numpy.copy(d_kplus1)
   simultaneousReduce(A, B)

   dimKChains = A.shape[1]
   kernelDim = dimKChains - numPivotCols(A)
   imageDim = numPivotRows(B)

   return kernelDim - imageDim

And there we have it! We’ve finally tackled the beast, and written a program to compute algebraic features of a topological space.

The reader may be curious as to why we didn’t come up with a more full-bodied representation of a simplicial complex and write an algorithm which accepts a simplicial complex and computes all of its homology groups. We’ll leave this direct approach as a (potentially long) exercise to the reader, because coming up in this series we are going to do one better. Instead of computing the homology groups of just one simplicial complex using by repeating one algorithm many times, we’re going to compute all the homology groups of a whole family of simplicial complexes in a single bound. This family of simplicial complexes will be constructed from a data set, and so, in grandiose words, we will compute the topological features of data.

If it sounds exciting, that’s because it is! We’ll be exploring a cutting-edge research field known as persistent homology, and we’ll see some of the applications of this theory to data analysis.

Until then!

Row Reduction Over A Field

We’re quite eager to get to applications of algebraic topology to things like machine learning (in particular, persistent homology). Even though there’s a massive amount of theory behind it (and we do plan to cover some of the theory), a lot of the actual computations boil down to working with matrices. Of course, this means we’re in the land of linear algebra; for a refresher on the terminology, see our primers on linear algebra.

In addition to applications of algebraic topology, our work with matrices in this post will allow us to solve important optimization problems, including linear programming. We will investigate these along the way.

Matrices Make the World Go Round

Fix two vector spaces V, W of finite dimensions m,n (resp), and fix bases for both spaces. Recall that an n \times m matrix uniquely represents a linear map V \to W, and moreover all such linear maps arise this way. Our goal is to write these linear maps in as simple a way as possible, so that we can glean useful information from them.

The data we want from a linear map f are: a basis for the kernel of f, a basis for the image of f, the dimensions of both, and the eigenvalues and eigenvectors of f (which is only defined if in addition V=W). As we have already seen, just computing the latter gives us a wealth of information about things like the internet and image similarity. In future posts, we will see how the former allows us to infer the shape (topology) of a large data set.

For this article, we will assume our vector spaces lie over the field \mathbb{R}, but later we will relax this assumption (indeed, they won’t even be vector spaces anymore!). Luckily for us, the particular choice of bases for V and W is irrelevant. The function f is fixed, and the only thing which changes is the matrix representation of f with respect to our choice of bases. If we pick the right bases, then the pieces of information above are quite easy to determine. Rigorously, we say two n \times m matrices A,B are row equivalent if there exists an invertible matrix M with B = MA. The reader should determine what the appropriate dimensions of this matrix are (hint: it is square). Moreover, we have the following proposition (stated without proof) which characterizes row equivalent matrices:

Proposition: Two matrices are row equivalent if and only if they represent the same linear map V \to W up to a choice of basis for W.

The “row” part of row equivalence is rather suggestive. Indeed, it turns out that our desired form B can be achieved by a manipulation of the rows of A, which is equivalently the left multiplication by a special matrix M.

Reduced Row Echelon Form, and Elementary Matrices

Before we go on, we should describe the convenient form we want to find. The historical name for it is reduced row echelon form, and it imposes a certain form on the rows of a matrix A:

  • its nonzero rows are above all of its zero rows,
  • the leftmost entry in each nonzero row is 1,
  • the leftmost entry in a row is strictly to the left of the leftmost entry in the next row, and
  • the leftmost entry in each nonzero row is the only nonzero entry in its column.

The reduced row echelon form is canonical, meaning it is uniquely determined for any matrix (this is not obvious, but for brevity we refer the reader to an external proof). For example, the following matrix \mathbb{R}^5 \to \mathbb{R}^3 has the given reduced row echelon form:

\begin{pmatrix} 1 & 2 & 3 & 4 & 5 \\ 1 & 2 & 4 & 5 & 6 \\ 0 & 0 & 2 & 1 & 0 \end{pmatrix} \to \begin{pmatrix} 1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & -1 \\ 0 & 0 & 0 & 1 & 2 \end{pmatrix}

The entries which contain a 1 and have zeros in the corresponding row and column are called pivots, and this name comes from their use in the algorithm below.

Calling this a “form” of the original matrix is again suggestive: a matrix in reduced row echelon form is just a representation with respect to a different basis (in fact, just a different basis for the codomain W). To prove this, we will show a matrix is row equivalent to its reduced row echelon form. But before we do that, we should verify that the reduced row echelon form actually gives us the information we want.

For the rightmost matrix above, and assuming we know the correct choice of basis for W is w_1, \dots, w_k, we can determine a basis for the image quite easily. Indeed, if the jth column contains a single 1 in the ith row, then the vector f(v_j) = w_i is in the image of f. Moreover, if we do this for each nonzero row (and because each nonzero row has a pivot) we obtain a basis for the whole image of f as a subset of the w_1 , \dots, w_k. Indeed, they are linearly independent as they form a basis of W, and they span the image of f because each basis vector w_i which was not found in the above way corresponds to a row of all zeros. In other words, it is clear from the entries of the reduced row echelon form matrix that no vectors f(v) expand as a linear combination of the unchosen w_i.

To put a matrix into reduced row echelon form, we allow ourselves three elementary row operations, and give an algorithm describing in what order to perform them. The operations are:

  • swap the positions of any two rows,
  • multiply any row by a nonzero constant, and
  • add a nonzero multiple of one row to another row.

Indeed, we should verify that these operations behave well. We can represent each one by the left multiplication by an appropriate matrix. Swapping the ith and jth rows corresponds to the identity matrix with the same rows swapped. Multiplying the ith row by a constant c corresponds to the identity matrix with c in the i,i entry. And adding c times row j to row i corresponds to the identity matrix with a c in the i,j entry. We call these matrices elementary matrices, and any sequence of elementary row operations corresponds to left multiplication by a product of elementary matrices. As we will see by our algorithm below, these row operations are enough to put any matrix (with entries in a field) into reduced row echelon form, hence proving that all matrices are row equivalent to some matrix in reduced row echelon form.

Before we describe this algorithm, we should make one important construction which will be useful in the future. Fixing the dimension n of our elementary matrices we note three things: the identity matrix is an elementary matrix, every elementary matrix is invertible, and the inverse of an elementary matrix is again an elementary matrix. In particular, every product of elementary matrices is invertible (the product of the inverses in reverse order), and so we can describe the group generated by all elementary matrices. We call this group the general linear group, denoted \textup{GL}_n( \mathbb{R}), and note that it has a very important place in the theory of Lie groups, which is (very roughly) the study of continuous symmetry. It has its name because we note that a matrix is in \textup{GL}_n(\mathbb{R}) if and only if its columns are linearly independent (and invertible). This happens precisely when it is row equivalent to the identity matrix. In other words, for any such matrix A there exists a product of elementary matrices E_1 \dots E_n such that E_1 \dots E_n A = I_n, and hence A = E_n^{-1} \dots E_1^{-1}. So we can phrase the question of whether a matrix is invertible as whether it is in \textup{GL}_n(\mathbb{R}), and answer it by finding its reduced row echelon form. So without further ado:

The Row Reduction Algorithm

Now that we’ve equipped ourselves with the right tools, let’s describe the algorithm which will transform any matrix into reduced row echelon form. We will do it straight in Python, explaining the steps along the way.

def rref(matrix):
   numRows = len(matrix)
   numCols = len(matrix[0])

   i,j = 0,0

The first part is straightforward: get the dimensions of the matrix and initialize i,j to 0. i,j represent the current row and column under inspection, respectively. We then start a loop:

   while True:
      if i >= numRows or j >= numCols:
         break

      if matrix[i][j] == 0:
         nonzeroRow = i
         while nonzeroRow < numRows and matrix[nonzeroRow][j] == 0:
            nonzeroRow += 1

      if nonzeroRow == numRows:
         j += 1
         continue

Here we check the base cases: if our indices have exceeded the bounds of our matrix, then we are done. Next, we need to find a pivot and put it in the right place, and essentially we work by induction on the columns. Since we are working over a field, any nonzero element can be a pivot, as we my just divide the entire row by the value of the leftmost entry to get a 1 in the right place. We just need to find a row with a nonzero value, and we prefer to pick the row which has the leftmost nonzero value, and if there are many rows with that property we pick the one in the row with the smallest index. In other words, we fix the leftmost column, try to find a pivot there by scanning downwards, and if we find none, we increment the column index and begin our search again. Once we find it, we may swap the two rows and save the pivot:

      temp = matrix[i]
      matrix[i] = matrix[nonzeroRow]
      matrix[nonzeroRow] = temp

      pivot = matrix[i][j]
      matrix[i] = [x / pivot for x in matrix[i]]

Once we have found a pivot, we simply need to eliminate the remaining entries in the column. We know this won’t affect any previously inspected columns, because by the inductive hypothesis any entries which are to the left of our pivot are zero.

      for otherRow in range(0, numRows):
         if otherRow == i:
         continue
         if matrix[otherRow][j] != 0:
            matrix[otherRow] = [y - matrix[otherRow][j]*x
             for (x,y) in zip(matrix[i], matrix[otherRow])]

      i += 1; j+= 1

   return matrix

After zeroing out the entries in the jth column, we may start to look for the next pivot, and since it can’t be in the same column or row, we may restrict our search to the sub-matrix starting at the i+1, j+1 entry. Once the while loop has terminated, we have processed all pivots, and we are done.

We encourage the reader to work out a few examples of this on small matrices, and modify our program to print out the matrix at each step of modification to verify the result. As usual, the reader may find the entire program on this blog’s Github page.

From here, determining the information we want is just a matter of reading the entries of the matrix and presenting it in the desired way. To determine the change of basis for W necessary to write the matrix as desired, one may modify the above algorithm to accept a second square matrix Q whose rows contain the starting basis (usually, the identity matrix for the standard basis vectors), and apply the same elementary row operations to Q as we do to A. The reader should try to prove that this does what it should, and we leave any further notes to a discussion in the comments.

Next time, we’ll relax the assumption that we’re working over a field. This will involve some discussion of rings and R-modules, but in the end we will work over very familiar number rings, like the integers \mathbb{Z} and the integers modulo n, and our work will be similar to the linear algebra we all know and love.

Until then!