Polyhedra, vertices and the fundamental theorem of LP

Lesson 2 showed in pictures that a two-variable linear program, if it has an optimum, has one at a corner of the feasible polygon. This lesson proves the statement in any number of dimensions. The feasible set of an LP is a polyhedron, and its corners can be described in three ways: geometrically as extreme points, through optimization as vertices, and algebraically as basic feasible solutions. Proving that the three descriptions agree is the bridge between the geometry of LP and the algebra of the simplex method. The payoff is the fundamental theorem of linear programming: a feasible LP is either unbounded or has an optimal basic feasible solution. Since there are only finitely many basic feasible solutions, solving an LP is in principle a finite search, and the remaining question, taken up in lesson 4, is how to search cleverly.

Polyhedra

A polyhedron is a set of the form \(P = \{x \in \mathbb{R}^n : Ax \le b\}\), the intersection of finitely many closed half-spaces \(a_i^Tx \le b_i\), where \(a_i^T\) is row \(i\) of the \(m \times n\) matrix \(A\). A bounded polyhedron is called a polytope. The feasible set of every LP is a polyhedron, since equations can be written as pairs of inequalities. Polyhedra are closed, as intersections of closed sets, and convex, by the argument of lesson 2.

For a point \(x \in P\), the constraint \(i\) is active (or tight, or binding) at \(x\) if \(a_i^Tx = b_i\). Write \(I(x)\) for the set of active constraints. Active constraints are the ones that hold \(x\) in place: at an interior point no constraint is active and \(x\) can move freely in all \(n\) directions; on a face of dimension \(k\) the active constraints leave \(k\) degrees of freedom; and at a corner they pin \(x\) down completely.

Three descriptions of a corner

Let \(P\) be a polyhedron and \(x \in P\). We consider three properties.

Extreme point. \(x\) is not a proper convex combination of two other points of \(P\): there are no \(u \ne w\) in \(P\) and \(\lambda \in (0, 1)\) with \(x = \lambda u + (1 - \lambda)w\). Geometrically, \(x\) does not lie strictly inside any segment in \(P\).

Vertex. There is a vector \(c\) such that \(x\) is the unique minimizer of \(c^Ty\) over \(P\): \(c^Tx < c^Ty\) for every \(y \in P\) with \(y \ne x\). Optimization-wise, some objective singles out \(x\).

Basic feasible solution. Among the constraints active at \(x\) there are \(n\) whose normal vectors \(a_i\) are linearly independent. Algebraically, the active constraints determine \(x\) as the unique solution of a square nonsingular linear system.

Theorem. For a point \(x\) of a nonempty polyhedron \(P\), the three properties are equivalent.

Proof. Vertex implies extreme point. Let \(x\) be the unique minimizer of \(c^Ty\) on \(P\) and suppose \(x = \lambda u + (1 - \lambda)w\) with \(u \ne w\) in \(P\) and \(0 < \lambda < 1\). Then \(u \ne x\) and \(w \ne x\) (if \(u = x\), then \((1 - \lambda)(w - x) = 0\) forces \(w = x = u\)). By uniqueness \(c^Tu > c^Tx\) and \(c^Tw > c^Tx\), so \(c^Tx = \lambda c^Tu + (1 - \lambda)c^Tw > c^Tx\), a contradiction.

Extreme point implies basic feasible solution. Suppose the active normals \(\{a_i : i \in I(x)\}\) span a subspace of dimension less than \(n\). Then there is a nonzero \(d\) orthogonal to all of them: \(a_i^Td = 0\) for \(i \in I(x)\). For the inactive constraints \(a_i^Tx < b_i\), so for \(\varepsilon > 0\) small enough, \(a_i^T(x \pm \varepsilon d) < b_i\) still holds; for the active ones \(a_i^T(x \pm \varepsilon d) = b_i\). Hence \(x + \varepsilon d\) and \(x - \varepsilon d\) both lie in \(P\), are distinct, and average to \(x\). So \(x\) is not extreme.

Basic feasible solution implies vertex. Let \(c = \sum_{i \in I(x)} a_i\). For any \(y \in P\), \(c^Ty = \sum_{i \in I(x)} a_i^Ty \le \sum_{i \in I(x)} b_i = c^Tx\). So \(x\) maximizes \(c^Ty\), and equivalently minimizes \((-c)^Ty\). Equality holds only if \(a_i^Ty = b_i\) for all \(i \in I(x)\); since \(n\) of these equations have independent normals, they have the unique solution \(y = x\). So \(x\) is the unique maximizer, and a vertex for the objective \(-c\). \(\square\)

Since the three notions agree, we use the words vertex, extreme point and basic feasible solution interchangeably. An immediate consequence is finiteness: each basic feasible solution is the unique solution of \(n\) of the \(m\) constraint equations, and there are at most \(\binom{m}{n}\) ways to choose them, so a polyhedron has at most \(\binom{m}{n}\) vertices. For the polytope in the figure, with seven constraints in three variables, that bound is 35, while the actual count is 8.

Basic solutions in standard form

The simplex method works with the standard form \(P = \{x \in \mathbb{R}^n : Ax = b,\ x \ge 0\}\), where now \(A\) is \(m \times n\) with \(m \le n\), and we assume the rows of \(A\) are linearly independent (redundant equations can be deleted). Here the constraints are the \(m\) equations, always active, and the \(n\) sign constraints \(x_j \ge 0\). A point has \(n\) independent active constraints exactly when, in addition to the \(m\) equations, enough sign constraints are tight. This leads to a purely algebraic recipe.

Choose \(m\) indices \(B = \{B(1), \ldots, B(m)\}\) such that the corresponding columns of \(A\) form an invertible \(m \times m\) matrix, also called \(B\), the basis matrix. Set the other \(n - m\) variables, the nonbasic ones, to zero, and solve \(Bx_B = b\) for the basic variables. The resulting \(x\) is a basic solution; if also \(x_B \ge 0\), it is a basic feasible solution.

Theorem. A point \(x \in P\) is a basic feasible solution if and only if the columns \(A_j\) with \(x_j > 0\) are linearly independent. In that case it arises from at least one basis \(B\) as above.

Proof. If \(x\) comes from a basis, the positive entries are among the basic ones, whose columns are independent. Conversely suppose the columns \(\{A_j : x_j > 0\}\) are independent. Since \(A\) has rank \(m\), we can extend them to \(m\) independent columns \(B\). The nonbasic entries of \(x\) are zero, so \(Bx_B = Ax = b\), and \(x\) is the basic solution for \(B\). The active constraints at \(x\) are the \(m\) equations and the \(n - m\) nonbasic sign constraints; written as rows of an \(n \times n\) matrix, after reordering columns this matrix is \(\begin{pmatrix} B & N \\ 0 & I \end{pmatrix}\), which is invertible because \(B\) is. So \(x\) is basic feasible in the earlier sense. \(\square\)

A basic feasible solution is degenerate if some basic variable is zero. Then more than \(n\) constraints are active, and the same vertex comes from several bases. Degeneracy is harmless geometrically but causes trouble for the simplex method, as lesson 5 explains.

The fundamental theorem

Theorem (fundamental theorem of linear programming). Consider minimizing \(c^Tx\) over \(P = \{x : Ax = b,\ x \ge 0\}\). (a) If \(P\) is nonempty, it has a basic feasible solution. (b) If the LP has a finite optimal value, then that value is attained at a basic feasible solution. (c) Otherwise, if \(P\) is nonempty, the optimal value is \(-\infty\).

Proof. The heart of the proof is a reduction step. Let \(x \in P\) and let \(S = \{j : x_j > 0\}\) be its support. If the columns \(A_j\), \(j \in S\), are independent, \(x\) is already a basic feasible solution. Otherwise there is \(d \ne 0\) with \(Ad = 0\) and \(d_j = 0\) for \(j \notin S\). Replacing \(d\) by \(-d\) if necessary we may assume \(c^Td \le 0\), and when \(c^Td = 0\) we may also assume some \(d_j < 0\).

Points \(x + td\) satisfy \(A(x + td) = b\) for every \(t\). If \(c^Td < 0\) and \(d \ge 0\), then \(x + td \ge 0\) for all \(t \ge 0\) and \(c^T(x + td) \to -\infty\): the LP is unbounded. Otherwise some \(d_j < 0\), necessarily with \(j \in S\). Let \(t^* = \min\{-x_j/d_j : d_j < 0\}\), which is positive because \(x_j > 0\) on \(S\). The point \(x' = x + t^*d\) is still nonnegative, satisfies \(c^Tx' \le c^Tx\), and has at least one more zero entry than \(x\), since the index achieving the minimum becomes zero while no entry outside \(S\) changes. Its support is strictly smaller.

Repeating the step, after at most \(n\) iterations we either discover that the LP is unbounded or reach a basic feasible solution \(\hat x\) with \(c^T\hat x \le c^Tx\). This proves (a). For (b) and (c), suppose the LP is not unbounded, so the reduction always ends at a basic feasible solution. There are finitely many basic feasible solutions; let \(x^*\) be the best of them. For every feasible \(x\), the reduction yields a basic feasible \(\hat x\) with \(c^Tx^* \le c^T\hat x \le c^Tx\). So \(x^*\) is optimal and the optimal value is finite and attained. If instead the reduction ever finds an unbounded direction, the optimal value is \(-\infty\). \(\square\)

The theorem has a hidden hypothesis that deserves comment. A general polyhedron need not have any vertex: a half-plane \(\{(x, y) : y \ge 0\}\) has none, and minimizing \(y\) over it has an optimal value of 0 attained along a whole line, with no corner at all. A nonempty polyhedron has a vertex exactly when it contains no line, and standard form polyhedra never contain a line, because \(x + td \ge 0\) for all real \(t\) forces \(d = 0\). This is one reason the standard form is convenient. Converting free variables to differences of nonnegative ones introduces corners where there were none.

A related classical result, due to Hermann Minkowski for polytopes and extended by Hermann Weyl in 1935, says that a polytope is the convex hull of its finitely many vertices, and conversely the convex hull of finitely many points is a polytope. The two descriptions, by inequalities and by vertices, can be wildly different in size: the \(n\)-dimensional cube has \(2n\) facets but \(2^n\) vertices.

Adjacent vertices and edges

Two basic feasible solutions are adjacent if their bases differ in exactly one index; geometrically, they are joined by an edge of \(P\), a segment along which \(n - 1\) independent constraints stay active. Moving from a vertex to an adjacent one means releasing one active constraint and following the edge until a new constraint becomes active. In the standard form language, one nonbasic variable increases from zero, the basic variables adjust to keep \(Ax = b\), and the move stops when some basic variable reaches zero. That is precisely a pivot of the simplex method. The polytope in the figure has 8 vertices and 12 edges, and each vertex has exactly three neighbors, as expected when no vertex is degenerate in three dimensions.

Enumerating bases

The program puts the shop model of lesson 2 in standard form with slacks \(s_1, s_2, s_3\), so \(m = 3\) and \(n = 5\). It tries all \(\binom{5}{3} = 10\) choices of three columns, solves for the basic variables in exact arithmetic, and sorts the results into dependent column sets, infeasible basic solutions and basic feasible solutions. It then runs the reduction from the proof of the fundamental theorem, starting from the interior point \((x, y) = (10, 10)\), and shows it arriving at a vertex with a better objective.

from fractions import Fraction as Fr
from itertools import combinations

# shop model in standard form: minimize c^T x, A x = b, x >= 0, variables x, y, s1, s2, s3
names = ["x", "y", "s1", "s2", "s3"]
A = [[2, 1, 1, 0, 0], [1, 1, 0, 1, 0], [1, 0, 0, 0, 1]]
b = [100, 80, 40]
c = [-30, -20, 0, 0, 0]
m, n = len(A), len(A[0])

def solve(M, rhs):
    """Gauss-Jordan elimination in exact arithmetic; None if M is singular."""
    k = len(M)
    T = [[Fr(v) for v in row] + [Fr(r)] for row, r in zip(M, rhs)]
    for col in range(k):
        piv = next((r for r in range(col, k) if T[r][col] != 0), None)
        if piv is None:
            return None
        T[col], T[piv] = T[piv], T[col]
        T[col] = [v / T[col][col] for v in T[col]]
        for r in range(k):
            if r != col and T[r][col] != 0:
                T[r] = [a - T[r][col] * p for a, p in zip(T[r], T[col])]
    return [row[-1] for row in T]

feasible = []
for basis in combinations(range(n), m):
    B = [[A[i][j] for j in basis] for i in range(m)]
    xB = solve(B, b)
    if xB is None:
        print(f"basis {[names[j] for j in basis]}: columns dependent, not a basis")
        continue
    x = [Fr(0)] * n
    for j, v in zip(basis, xB):
        x[j] = v
    ok = all(v >= 0 for v in x)
    z = sum(ci * xi for ci, xi in zip(c, x))
    print(f"basis {[names[j] for j in basis]}: (x, y) = ({x[0]}, {x[1]}), "
          f"slacks {[str(v) for v in x[2:]]}, {'feasible, value ' + str(z) if ok else 'infeasible'}")
    if ok:
        feasible.append((z, x))
print(f"{len(feasible)} basic feasible solutions; best value {min(feasible)[0]} at "
      f"(x, y) = ({min(feasible)[1][0]}, {min(feasible)[1][1]})")

# the reduction in the proof of the fundamental theorem: start inside, end at a vertex
def null_vector(cols):
    """A nonzero d with A d = 0 supported on cols, or None if those columns are independent."""
    T = [[Fr(A[i][j]) for j in cols] for i in range(m)]
    pivots, r = [], 0
    for k in range(len(cols)):
        p = next((i for i in range(r, m) if T[i][k] != 0), None)
        if p is None:
            continue
        T[r], T[p] = T[p], T[r]
        T[r] = [v / T[r][k] for v in T[r]]
        for i in range(m):
            if i != r and T[i][k] != 0:
                T[i] = [a - T[i][k] * q for a, q in zip(T[i], T[r])]
        pivots.append(k); r += 1
    free = [k for k in range(len(cols)) if k not in pivots]
    if not free:
        return None
    f = free[0]
    d = [Fr(0)] * n
    d[cols[f]] = Fr(1)
    for row, k in enumerate(pivots):
        d[cols[k]] = -T[row][f]
    return d

x = [Fr(10), Fr(10), Fr(70), Fr(60), Fr(30)]       # feasible, all five entries positive
value = lambda x: sum(ci * xi for ci, xi in zip(c, x))
print("start", [str(v) for v in x], " value", value(x))
while (d := null_vector([j for j in range(n) if x[j] > 0])) is not None:
    cd = sum(ci * di for ci, di in zip(c, d))
    if cd > 0 or (cd == 0 and all(v >= 0 for v in d)):
        d = [-v for v in d]                   # objective must not increase; need a negative entry
    if all(v >= 0 for v in d):
        print("objective decreases forever along d: unbounded"); break
    t = min(-x[j] / d[j] for j in range(n) if d[j] < 0)   # step until an entry hits zero
    x = [xi + t * di for xi, di in zip(x, d)]
    print("moved to", [str(v) for v in x], " value", value(x))
print("reached a basic feasible solution with support", [names[j] for j in range(n) if x[j] > 0])

basis ['x', 'y', 's1']: (x, y) = (40, 40), slacks ['-20', '0', '0'], infeasible
basis ['x', 'y', 's2']: (x, y) = (40, 20), slacks ['0', '20', '0'], feasible, value -1600
basis ['x', 'y', 's3']: (x, y) = (20, 60), slacks ['0', '0', '20'], feasible, value -1800
basis ['x', 's1', 's2']: (x, y) = (40, 0), slacks ['20', '40', '0'], feasible, value -1200
basis ['x', 's1', 's3']: (x, y) = (80, 0), slacks ['-60', '0', '-40'], infeasible
basis ['x', 's2', 's3']: (x, y) = (50, 0), slacks ['0', '30', '-10'], infeasible
basis ['y', 's1', 's2']: columns dependent, not a basis
basis ['y', 's1', 's3']: (x, y) = (0, 80), slacks ['20', '0', '40'], feasible, value -1600
basis ['y', 's2', 's3']: (x, y) = (0, 100), slacks ['0', '-20', '40'], infeasible
basis ['s1', 's2', 's3']: (x, y) = (0, 0), slacks ['100', '80', '40'], feasible, value 0
5 basic feasible solutions; best value -1800 at (x, y) = (20, 60)
start ['10', '10', '70', '60', '30']  value -500
moved to ['10', '70', '10', '0', '30']  value -1700
moved to ['20', '60', '0', '0', '20']  value -1800
reached a basic feasible solution with support ['x', 'y', 's3']

Of the ten column triples, one is dependent (the columns of \(y\), \(s_1\) and \(s_2\) all have a zero in the demand row, so they span only a plane), four give basic solutions with a negative entry, and five give the five corners of the pentagon, matching the picture in lesson 2. The reduction from the interior point takes two steps: first it fills the finishing hours, then the carpentry hours, arriving at the optimal vertex \((20, 60)\). In general the reduction only guarantees some vertex no worse than the start, not the optimal one; that it found the optimum here is luck. The simplex method adds the missing ingredient, a test that tells us when a vertex is optimal and which neighbor to move to when it is not.