Algorithm Design and Analysis

Polytopes and Integrality

17 of 17 · 140 minLabLP Duality

Polytopes and Integrality

The question

A small agency has three workers, A, B and C, and three jobs, P, Q and R. Each worker takes at most one job and each job needs at most one worker. The value of each allowed assignment is known:

5 8 7 6 3 4 6 A B C P Q R
The hook: workers A, B, C on the left, jobs P, Q, R on the right. A line is an allowed assignment, labelled with its value.

The best assignment is easy to find by hand. The interesting part is what happens when you hand the problem to a linear-programming solver. Give every allowed pair a variable xe, meaning "how much of worker u goes to job v", and ask for

max∑ewexes.t.∑e∋vxe≤1  for every worker and job v,x≥0.

Nothing in this LP says that xe must be 0 or 1. Will the solver ever hand back half a worker? For this problem, never. For a problem that looks almost the same, it will. This module explains both answers, and on the way what the simplex method does and how long it can take.

You need module 15's LPs in the form maxc·x, Ax≤b, x≥0 and its tableau simplex, and module 11's matchings. Everything else is built here.

Rounding, and enumerating

Two obvious ways to get whole numbers out of a fractional solver both waste something.

Round the answer. Take three people u,v,w, any two of whom may be paired, each pair worth 1. The LP's best point turns out to be 12 on all three pairs, with value 3/2. Rounding every 12 up puts every person in two pairs, which is infeasible. Rounding down gives nothing, value 0, while one real pair is worth 1.

Enumerate. List every matching, or branch on every fractional variable (branch and bound). Three workers and three jobs have at most 3!=6 full assignments, but n and n have n!, and branching can double the work with each variable.

from fractions import Fraction as F
import itertools, math

HOOK = [("A", "P", 5), ("A", "Q", 8), ("B", "P", 7), ("B", "Q", 6),
        ("B", "R", 3), ("C", "Q", 4), ("C", "R", 6)]

def best_matching(edges):
    """The oracle: skip the first edge, or take it and drop every edge touching it."""
    if not edges:
        return 0
    (u, v, w), rest = edges[0], edges[1:]
    return max(best_matching(rest),
               w + best_matching([e for e in rest if not {u, v} & {e[0], e[1]}]))

assert best_matching(HOOK) == 21            # A-Q 8, B-P 7, C-R 6
assert math.factorial(10) == 3628800        # full assignments of 10 workers to 10 jobs

The waste is that neither approach uses what the LP already knows. For the hook, the solver's answer is already whole, and the reason is geometric.

Corners are basic feasible solutions

The feasible region P={x:Ax≤b, x≥0} is an intersection of half-spaces, a convex polyhedron; when it is bounded it is a polytope. A feasible point x is a vertex (a corner) if it is not the midpoint of two different feasible points. Module 15 used, without proof, that an LP with an optimum has an optimal vertex. Here is the algebraic description of vertices that makes that usable.

Call a constraint tight at x if it holds with equality: a row with ai·x=bi, or a bound with xj=0.

Theorem (vertices are basic feasible solutions)

A feasible point of P⊆Rn is a vertex if and only if the constraints tight at it include n linearly independent ones.

Proof. If the tight constraints have rank less than n, some direction d≠0 is orthogonal to all of them. Moving to x±εd keeps every tight constraint tight, and for a small enough ε>0 it keeps every slack constraint slack, so both points are feasible and x is their midpoint: not a vertex. Conversely, suppose x=(y+z)/2 with y≠z feasible. If a·x=β is tight, then a·y≤β and a·z≤β average to β, so both are equalities. So d=y−z≠0 satisfies a·d=0 for every tight constraint, and their rank is less than n. ◻

After adding a slack variable to each row, "rank n among the tight constraints" says the same thing as "the non-zero variables have linearly independent columns", which is module 15's basic feasible solution: the point a simplex tableau stands on. Three consequences follow. A vertex is the unique solution of n tight constraints, so by Cramer's rule its coordinates are ratios of determinants of submatrices of [A b]. There are at most (m+nn) vertices, one per choice of n tight constraints out of m+n. And that bound is the naive algorithm: try every choice.

Take the LP we will call L: maximize x+3y subject to x≤4, x+y≤6, −x+y≤3, y≤4 and x,y≥0. It has m+n=6 constraints, so (62)=15 pairs of lines to try:

ROWS = [(1, 0, 4), (1, 1, 6), (-1, 1, 3), (0, 1, 4), (-1, 0, 0), (0, -1, 0)]  # a*x + b*y <= r

def corners(rows):
    """Every intersection of two constraint lines that satisfies all the constraints."""
    solved, found = 0, set()
    for (a, b, r), (c, d, s) in itertools.combinations(rows, 2):
        det = a * d - b * c
        if det == 0:
            continue                          # parallel lines never meet
        solved += 1
        x, y = F(r * d - b * s, det), F(a * s - c * r, det)      # Cramer's rule
        if all(p * x + q * y <= t for p, q, t in rows):
            found.add((x, y))
    return solved, sorted(found)

solved, L_corners = corners(ROWS)
assert math.comb(6, 2) == 15 and solved == 13
assert L_corners == [(0, 0), (0, 3), (1, 4), (2, 4), (4, 0), (4, 2)]
assert max(x + 3 * y for x, y in L_corners) == 14          # at (2, 4)

Thirteen of the 15 pairs meet (two pairs are parallel), and only 6 of those 13 points are feasible. In n dimensions the number of candidate systems grows like (m+nn), which is exponential in n. The simplex method looks at far fewer.

The simplex method walks from corner to corner

A simplex tableau always represents a basic feasible solution, so by the theorem the method only ever stands on vertices. A pivot swaps one tight constraint for another: the entering column leaves its bound, and the leaving row becomes tight. Geometrically, it slides along an edge of the polytope to an adjacent vertex, and the objective never goes down. At a degenerate vertex, where more than n constraints are tight, a pivot can swap constraints without moving at all.

The pivot rule chooses the entering column among those with a negative reduced cost. Bland's rule takes the smallest index; it never cycles, which is why the course's toolkit uses it. Dantzig's rule takes the most negative reduced cost, the steepest improvement per unit. Here is a minimal tableau simplex that records every vertex it stands on (for b≥0, so the origin is a starting vertex, and for bounded LPs):

def walk(c, A, b, rule="bland"):
    """Tableau simplex for max c.x, Ax <= b, x >= 0 with b >= 0. Returns the vertices visited."""
    m, n = len(A), len(c)
    T = [[F(v) for v in A[i]] + [F(int(i == k)) for k in range(m)] + [F(b[i])] for i in range(m)]
    T.append([F(-v) for v in c] + [F(0)] * (m + 1))           # the objective row
    B = [n + i for i in range(m)]                               # the slacks start basic
    def point():
        x = [F(0)] * (n + m)
        for i, j in enumerate(B):
            x[j] = T[i][-1]
        return tuple(x[:n])
    path = [point()]
    while True:
        neg = [j for j in range(n + m) if T[-1][j] < 0]
        if not neg:
            return path
        s = neg[0] if rule == "bland" else min(neg, key=lambda j: (T[-1][j], j))
        r = min((i for i in range(m) if T[i][s] > 0), key=lambda i: (T[i][-1] / T[i][s], B[i]))
        T[r] = [v / T[r][s] for v in T[r]]
        for i in range(m + 1):
            if i != r:
                T[i] = [a - T[i][s] * p for a, p in zip(T[i], T[r])]
        B[r] = s
        path.append(point())

cL, AL, bL = [1, 3], [[1, 0], [1, 1], [-1, 1], [0, 1]], [4, 6, 3, 4]
assert walk(cL, AL, bL) == [(0, 0), (4, 0), (4, 2), (2, 4)]
assert walk(cL, AL, bL, "dantzig") == [(0, 0), (0, 3), (1, 4), (2, 4)]
assert all(p in L_corners for p in walk(cL, AL, bL) + walk(cL, AL, bL, "dantzig"))
1 2 3 4 5 1 2 3 4 5 x + y = 6 −x + y = 3 x = 4 y = 4 (0,0) 0 (4,0) 4 (4,2) 10 (2,4) 14 (1,4) 13 (0,3) 9
The region L: each corner is labelled with its coordinates and its objective x + 3y. Purple: Bland's walk along the bottom and the right, (0,0) → (4,0) → (4,2) → (2,4). Teal: Dantzig's walk up the left, (0,0) → (0,3) → (1,4) → (2,4). Three pivots each, and neither ever leaves the boundary.

Both walks visit 4 of L's 6 corners. Bland's rule enters x first because it has the smaller index, and so it goes the long way along the bottom. Dantzig's rule enters y, whose coefficient 3 is larger, and climbs the left side. On L the two cost the same.

Predict: Dantzig's rule takes the steepest edge at every step. Can it ever take more pivots than a rule that ignores slopes?

Yes. Steepest per unit of the entering variable says nothing about how long the edge is or where it leads. The next section builds a family on which Dantzig's rule visits every single corner.

How long can the walk be?

Klee and Minty (1972) squashed an n-dimensional cube so that the objective increases along a path through all 2n corners, and Dantzig's rule follows that path. One version is

max∑j<n2n−1−jxjs.t.∑j<i2i−j+1xj+xi≤5i+1  (i=0,…,n−1),x≥0.

In two dimensions it is max2x0+x1 with x0≤5 and 4x0+x1≤25:

5 5 10 15 20 25 (0,0) 0 (5,0) 10 (5,5) 15 (0,25) 25
The Klee–Minty square (n = 2): a square squashed so that the objective 2x₀ + x₁ rises around all four corners. Dantzig's rule starts along x₀, the steeper direction, and visits every corner: 3 pivots, optimum 25. The grey lines are x₀ = 5 and 4x₀ + x₁ = 25.
def klee_minty(n):
    c = [2 ** (n - 1 - j) for j in range(n)]
    A = [[2 ** (i - j + 1) if j < i else int(i == j) for j in range(n)] for i in range(n)]
    return c, A, [5 ** (i + 1) for i in range(n)]

assert walk(*klee_minty(2), "dantzig") == [(0, 0), (5, 0), (5, 5), (0, 25)]
path3 = walk(*klee_minty(3), "dantzig")
assert path3 == [(0, 0, 0), (5, 0, 0), (5, 5, 0), (0, 25, 0),
                 (0, 25, 25), (5, 5, 65), (5, 0, 85), (0, 0, 125)]
assert [4 * a + 2 * b + c for a, b, c in path3] == [0, 20, 30, 50, 75, 95, 105, 125]

In three dimensions the walk visits all 8 corners of the squashed cube, and the objective rises at every step: 0, 20, 30, 50, 75, 95, 105, 125. The pattern is recursive. The walk solves the (n−1)-dimensional cube on the "floor" xn−1=0, steps up once, then walks the floor's cube backwards on the "ceiling". So P(n)=2P(n−1)+1 pivots, which gives 2n−1. That recursion is Klee and Minty's theorem for their construction. For this particular family the course verifies the count by running it, as the next section shows.

Other rules do not escape. Exponential examples are known for Bland's rule and many other common pivot rules, and whether any pivot rule needs only polynomially many pivots is open.

LPs are still solvable in polynomial time, by methods that do not walk along edges. The ellipsoid method (Khachiyan, 1979) and interior-point methods (Karmarkar, 1984) run in time polynomial in the number of bits needed to write the input down. Both are stated here without proof. An interior-point method moves through the inside of the polytope and may end at a non-vertex point of the optimal face; a final "crossover" step moves to a vertex. In practice simplex is usually fast. Spielman and Teng (2004) proved that its expected number of pivots is polynomial once the input is perturbed slightly at random (smoothed analysis). Production solvers offer both kinds of method.

Measure the claim

Pivots of the walk above on the Klee–Minty family, under both rules:

n 1 2 3 4 5 6 7 8
Dantzig 1 3 7 15 31 63 127 255
2n−1 1 3 7 15 31 63 127 255
Bland 1 3 5 9 15 25 41 67
dantzig = [len(walk(*klee_minty(n), "dantzig")) - 1 for n in range(1, 9)]
bland = [len(walk(*klee_minty(n), "bland")) - 1 for n in range(1, 9)]
assert dantzig == [2 ** n - 1 for n in range(1, 9)]
assert bland == [1, 3, 5, 9, 15, 25, 41, 67]
assert all(bland[k] == bland[k - 1] + bland[k - 2] + 1 for k in range(2, 8))
assert all(walk(*klee_minty(n), r)[-1][-1] == 5 ** n for n in range(1, 9) for r in ("bland", "dantzig"))

Dantzig's count is exactly 2n−1 on every cube measured. Bland's count is smaller, but each entry is the sum of the previous two plus one, so on these instances it grows by a factor of about 1.6 per dimension: still exponential, just with a smaller base. These are measurements on one family, n≤8. They illustrate the theorems but do not prove them.

When the corners are whole numbers

Back to the question. Run the same walk on the hook's matching LP:

def matching_lp(vertices, edges):
    A = [[int(v in (e[0], e[1])) for e in edges] for v in vertices]
    return [e[2] for e in edges], A, [1] * len(vertices)

hook_path = walk(*matching_lp(["A", "B", "C", "P", "Q", "R"], HOOK))
assert len(hook_path) - 1 == 9
assert all(v in (0, 1) for point in hook_path for v in point)     # every corner visited is 0/1
assert hook_path[-1] == (0, 1, 1, 0, 0, 0, 1)                      # A-Q, B-P, C-R
assert sum(e[2] * v for e, v in zip(HOOK, hook_path[-1])) == 21 == best_matching(HOOK)

Nine pivots, and every corner it stands on, not just the last one, is a 0/1 vector: a matching. That is no accident of this instance.

Theorem (bipartite matching polytope). If the graph is bipartite, every vertex of {x≥0:∑e∋vxe≤1 for all v} is a 0/1 vector.

The proof shows that any feasible point with a fractional edge is a midpoint, so it cannot be a vertex. Take the point with 12 on A–P, A–Q, B–P, B–Q and 1 on C–R. It is feasible, and its weight is (5+8+7+6)/2+6=19. Its fractional edges form the cycle A–P–B–Q–A. Put d=+1 on A–P and B–Q and d=−1 on A–Q and B–P, alternating around the cycle. Every vertex of the cycle meets one +1 edge and one −1 edge, so its load does not change under x±εd:

X = [F(1, 2), F(1, 2), F(1, 2), F(1, 2), 0, 0, 1]     # edge order as in HOOK
d = [1, -1, -1, 1, 0, 0, 0]                            # alternate around A-P-B-Q-A

def load(edges, x):
    total = {}
    for (u, v, _), xe in zip(edges, x):
        total[u] = total.get(u, 0) + xe
        total[v] = total.get(v, 0) + xe
    return total

def weight(edges, x):
    return sum(e[2] * xe for e, xe in zip(edges, x))

up = [a + F(1, 2) * b for a, b in zip(X, d)]
down = [a - F(1, 2) * b for a, b in zip(X, d)]
assert load(HOOK, X) == load(HOOK, up) == load(HOOK, down)       # every load unchanged
assert (up, down) == ([1, 0, 0, 1, 0, 0, 1], [0, 1, 1, 0, 0, 0, 1])
assert (weight(HOOK, X), weight(HOOK, up), weight(HOOK, down)) == (19, 17, 21)
Splitting along the cycle A–P–B–Q1x: weight 192x + d/2: weight 173x − d/2: weight 21
Top: the point, with its four fractional edges highlighted. Middle and bottom: the two points reached by pushing ½ around the cycle each way (edges at 1 highlighted). Both are matchings, and 19 = (17 + 21)/2.

Invariant

In a bipartite graph, a feasible point with a fractional edge lies in the middle of an even cycle, or of a path, of fractional edges, and pushing ±ε alternately along it keeps every load within its limit. So the point is a midpoint, not a vertex.

Proof in full. Let x be feasible with at least one fractional edge, and let F be the set of fractional edges.

  • F contains a cycle. The graph is bipartite, so the cycle is even, and signs +1,−1,+1,… around it close up consistently. Every vertex on the cycle meets one edge of each sign, so its load is unchanged. Take ε=minemin(xe,1−xe)>0 over the cycle's edges. Then both x±εd stay in [0,1] on every edge and keep every load, so both are feasible.
  • F is a forest. Take a longest path in it. Its interior vertices are balanced as before. An end vertex meets exactly one fractional edge, since otherwise the path could be extended (or would close a cycle). Its other edges are integral, and none can be 1, because together with the fractional edge its load would exceed 1. So its load equals that one edge's value, which is below 1, and moving that edge by ±ε with ε=minemin(xe,1−xe) keeps the load in [0,1].

Either way x=12(x+εd)+12(x−εd) with the two points different, so x is not a vertex. Every vertex is therefore 0/1, and a 0/1 feasible point is exactly a matching. ◻

The proof is also an algorithm. Push in whichever direction does not lower the weight, and go as far as feasibility allows, until some edge reaches 0 or 1. Each push makes at least one more edge whole and never makes a whole edge fractional. So at most |E| pushes turn any feasible point into a matching at least as heavy.

Where it breaks: the triangle

Now the three people u,v,w, any two of whom may pair up, each pair worth 1.

TRI = [("u", "v", 1), ("v", "w", 1), ("u", "w", 1)]
c3, A3, b3 = matching_lp(["u", "v", "w"], TRI)
tri_path = walk(c3, A3, b3)
assert tri_path[-1] == (F(1, 2), F(1, 2), F(1, 2)) and sum(tri_path[-1]) == F(3, 2)
assert tri_path[1] == tri_path[2]              # one degenerate pivot: it swapped, did not move
assert best_matching(TRI) == 1

def det3(M):
    (a, b, c), (d, e, f), (g, h, i) = M
    return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)

assert A3 == [[1, 0, 1], [1, 1, 0], [0, 1, 1]] and det3(A3) == 2
1/2 1/2 1/2 u v w
The triangle's LP optimum: ½ on every edge, value 3/2. Every vertex carries load exactly 1, yet no matching has more than one edge.

The point (12,12,12) is a genuine vertex. All three vertex constraints are tight there, and their matrix has determinant 2≠0, so they are independent: rank 3 in 3 dimensions. The simplex method stops on it with value 3/2, while the best real pairing is worth 1. The walk also makes one degenerate pivot, changing the basis without moving.

The proof fails at exactly one step. The fractional edges form a triangle, an odd cycle. Alternating +,−,+ around it puts two + edges on one vertex, whose load rises by 2ε. There is no direction to split along. Cramer's rule shows where the 12 comes from: a vertex's coordinates are ratios with the determinant, here 2, in the denominator.

Integrality for every graph needs more constraints. Edmonds (1965) showed that adding ∑e⊆Sxe≤(|S|−1)/2 for every odd set S of vertices makes the polytope integral for all graphs. That theorem is stated here, not proved.

Predict: add a fourth person z who can pair only with w. Is the LP still fractional?

No. Now xwz=1 and xuv=1 is a matching of weight 2, and the LP cannot beat 2: each of the four people has load at most 1 and every pair uses two people, so the value is at most 4/2=2. An odd cycle only causes a gap when the LP can profit from it.

Total unimodularity

The determinant is the whole story, and it has a name. A matrix is totally unimodular (TU) if every square submatrix has determinant −1, 0 or 1.

Theorem. If A is TU and b is integral, every vertex of {x:Ax≤b, x≥0} is integral. Proof. A vertex is the unique solution of n independent tight constraints: a non-singular square system Mx=b′ whose rows come from A or from the identity (the tight bounds xj=0). Expanding along the identity rows shows detM is ± a square subdeterminant of A, so detM=±1. Cramer's rule divides integers by ±1. ◻

Lemma. The vertex–edge incidence matrix of a bipartite graph is TU. Proof sketch. Induct on the size of a square submatrix M. If some column of M is zero, detM=0. If some column has a single 1, expand along it to a smaller submatrix, which is TU by induction. Otherwise every column has two 1s, one in a row of each side. Then the rows of one side and the rows of the other side both sum to the all-ones vector, so the rows are dependent and detM=0. ◻

The same argument works for directed incidence matrices, with one +1 and one −1 per column, which is why max flow and min-cost flow with integer capacities have integral optimal solutions (modules 11–13). The triangle's matrix is not TU, since its determinant is 2. Neither is a general coefficient matrix such as (2112), whose determinant is 3.

Complexity achieved. With a TU matrix and integral data, a yes-or-no problem is solved as an LP with no rounding step: any method that returns an optimal vertex returns an integral answer. A polynomial LP method plus a final move to a vertex gives a polynomial algorithm for weighted bipartite matching, instead of n! enumeration. Simplex itself has no polynomial pivot guarantee (Klee–Minty).

A problem that looks different

A city's streets run north–south and east–west. A camera placed at an intersection watches both streets through it. The city wants every street watched by the fewest cameras, and it plans to ask a solver that allows fractional cameras. Can it trust the solver to answer in whole cameras? Nothing here mentions workers or jobs. The lab's last problem is a different one.

Practise

The lab starts from the algebra. You will write the test "is this point a vertex?" and watch it hold at every point the simplex method visits on a new LP. Then you predict and carry out one split of a fractional point along an alternating cycle, and turn the split into a procedure that makes any fractional point of a bipartite graph whole. After that you measure pivot counts on a second Klee–Minty family under two rules. The last problem does not say what it is.

Recap

You can now: decide whether a point is a vertex by the rank of its tight constraints; follow the simplex method as a walk along edges between vertices; prove that bipartite matching LPs have 0/1 vertices, and show that a triangle breaks this; state total unimodularity and use it to trust an LP's answer to be whole.

Invariant: in a bipartite graph, a fractional feasible point sits in the middle of an even alternating cycle or path of fractional edges, so it is not a vertex.

Complexity achieved: an integral optimum from one LP solve, with no rounding and no branching, against n! enumeration. Simplex can take 2n−1 pivots (Dantzig's rule on Klee–Minty, and exponentially many for Bland's rule on the measured cubes), while the ellipsoid and interior-point methods are polynomial in the input's bit length.

Failure mode: believing that "LP optima are vertices" means "LP optima are integral". The triangle's (12,12,12) is a vertex too. Integrality comes from the constraint matrix, and rounding a fractional answer can be infeasible or far from optimal.

In real software: the HiGHS solver, which SciPy's scipy.optimize.linprog calls, offers both a dual simplex method (method="highs-ds") and an interior-point method (method="highs-ipm"). Mixed-integer solvers such as HiGHS and SCIP solve LP relaxations inside branch and bound, so an integral relaxation ends the search at its first node.

Retrieval: module 11 proved that a network with whole-number capacities has a whole-number maximum flow without any LP. What was that argument, and which lemma above gives the same conclusion from the constraint matrix?

Check yourself

  1. Why is the point with 12 on the cycle A–P–B–Q not a vertex, and exactly where does the argument use that the graph is bipartite?
  2. Run the hook's matching LP with Dantzig's rule instead of Bland's. Could the walk ever stand on a point with a fractional coordinate? Why, or why not?
  3. On the Klee–Minty cube with n=8, compare Dantzig's rule with Bland's rule. What do the counts show, and what do they not show?

Topic list modeled on the public schedule of CMU 15-451/651 Algorithm Design and Analysis (Fall 2025). All lecture text, proofs, examples, code, and exercises are original to SciMigo; no CMU course materials are reproduced. This course is independent and not affiliated with or endorsed by Carnegie Mellon University.

Your turn

Now try it yourself

Write the code, press Run, and scrub through the frames your program draws. The checks tell you when it works.

Loading the lab…