Linear Programming

The question

A small workshop makes benches and shelves. A bench earns 7 and a shelf earns 5. A bench takes 2 hours of cutting and 1 hour of finishing, and a shelf takes 1 hour of cutting and 2 of finishing. This week there are 14 cutting hours and 16 finishing hours, and the market will take at most 6 benches. Writing $x$ for benches and $y$ for shelves, the owner wants

$$max7x+5y\text{subject to}2x+y\le 14,x+2y\le 16,x\le 6,x,y\ge 0.$$

Every rule is a sum of the unknowns times fixed numbers, compared with a fixed number, and so is the profit. A problem of that shape is a linear program (LP). This module shows why its answer sits at a corner, how the simplex method walks from corner to corner, and how a diet, a knapsack, a maximum flow and a shortest path turn out to be LPs too.

What is the best plan, and how would the owner convince her partner that no plan earns more? This module answers the first question; module 16 answers the second.

You need two things from earlier modules: flows and cuts (module 11) for one of the examples, and the idea from module 14 that a game's value is the best of the worst cases. Nothing else.

Trying plans, and trying corners

The obvious approach is to try plans. With whole numbers there are only 7 × 9 candidates here, and the best whole-number plan is 4 benches and 6 shelves, profit 58. But the unknowns of an LP are real numbers. Change the finishing hours from 16 to 17 and the best plan becomes $x=11/3$, $y=20/3$, profit 59. No grid of whole numbers finds it, and no grid, however fine, proves that nothing between its points does better.

from fractions import Fraction as F
from math import comb

def best_whole_plan(finishing):
    return max((7 * x + 5 * y, x, y) for x in range(7) for y in range(9)
               if 2 * x + y <= 14 and x + 2 * y <= finishing)

assert best_whole_plan(16) == (58, 4, 6)
assert best_whole_plan(17) == (58, 4, 6)   # the grid is stuck at 58
x, y = F(11, 3), F(20, 3)                  # but this plan fits
assert 2 * x + y == 14 and x + 2 * y == 17 and x <= 6
assert 7 * x + 5 * y == 59

The better naive approach uses a fact shown below: some corner of the region is optimal. A corner of a region in two variables is where two of its boundary lines meet. The workshop has five boundary lines (three rules and $x=0$, $y=0$), so there are at most $\left(\genfrac{}{}{0ex}{}{5}{2}\right)=10$ pairs to intersect. Solve each pair, throw away the points that break a rule, and keep the best of the rest. In general, with $n$ variables and $m$ rules, a corner is where $n$ of the $m+n$ boundaries are tight, so this vertex enumeration solves up to $\left(\genfrac{}{}{0ex}{}{m+n}{n}\right)$ linear systems. For six variables and six rules that is 924 systems; for 30 and 30 it is about ${10}^{17}$. The waste: most systems land on a point that breaks a rule, or on a corner far from the best one.

The model

A linear program has real unknowns $x\in {\mathbb{R}}^n$, one linear objective $c·x$ to maximize or minimize, and finitely many linear constraints, each $a·x\le \beta$, $a·x\ge \beta$ or $a·x=\beta$. The points satisfying every constraint form the feasible region. Each constraint is a half-space (or a hyperplane), and an intersection of half-spaces is convex: the segment between two feasible points is feasible. Exactly one of three things is true of an LP:

• it is infeasible: the region is empty;
• it is unbounded: the objective grows without limit on the region;
• it has an optimum, a feasible point that no other feasible point beats.

This theorem is stated here; the simplex method below is a constructive proof of it. We will write every LP in one form, $maxc·x$ subject to $Ax\le b$ and $x\ge 0$, and do all arithmetic in exact fractions, so an answer of $59$ or $91/6$ is exact and never $58.99999$. What we count is pivots (defined below) and, for the naive approach, linear systems solved.

Slide the objective line

In two variables you can see the answer. The points of profit $z$ lie on the line $7x+5y=z$, and changing $z$ slides that line without turning it. Start at $z=0$ (through the origin) and slide it in the direction of $(7,5)$. The last value at which it still touches the region is the optimum. On the workshop that happens at $(4,6)$ with $z=58$: the dashed line in the figure touches the region only there.

Sketch. Take an optimal point $x$ that is not a vertex. Then there is a direction $d\ne 0$ with $x+d$ and $x-d$ both feasible (a point on an edge or inside the region can move both ways). The two values $c·(x\pm d)$ are at most $c·x$ and average to it, so both equal it: $c·d=0$. Move from $x$ along $d$ (or $-d$) until a constraint that was loose becomes tight. The value does not change, and one more constraint is tight. After at most $n$ such moves, $n$ independent constraints are tight, which is a vertex. The step that needs care is that the move stops at all; the full proof, with the precise definition of a vertex, is in module 17. $◻$

The corners are the same whatever the prices: only which corner the line touches last depends on them. At $(4,6)$ the cutting and finishing lines are tight, with slopes $-2$ and $-1/2$. The profit line has slope $-7/5$, between the two, which is why it touches that corner last.

CORNERS = [(0, 0), (6, 0), (6, 2), (4, 6), (0, 8)]
RULES = [([2, 1], 14), ([1, 2], 16), ([1, 0], 6)]

def slacks(x, y):
    """How much of each rule is left over; 0 means the rule is tight."""
    return [r - (a * x + b * y) for (a, b), r in RULES] + [x, y]

for x, y in CORNERS:
    s = slacks(x, y)
    assert min(s) >= 0 and s.count(0) >= 2   # feasible, two tight
assert [7 * x + 5 * y for x, y in CORNERS] == [0, 42, 52, 58, 40]
assert max((7 * x + 15 * y, (x, y)) for x, y in CORNERS) == (120, (0, 8))
assert 7 * 4 + 14 * 6 == 14 * 8 == 112     # tie at shelf profit 14

Walking the corners: the simplex method

Vertex enumeration visits every corner, most of them useless. The simplex method starts at one corner and repeatedly moves to a neighbouring corner that is at least as good, until no neighbour is better. In two variables the neighbours of a corner are the two corners at the ends of its two edges.

To do this algebraically, give each rule a slack variable: $s_1=14-2x-y$, $s_2=16-x-2y$, $s_3=6-x$, all required to be $\ge 0$. Number the variables $x,y,s_1,s_2,s_3$ as $0,1,2,3,4$. At a corner, two of the five are zero (the tight boundaries); the other three are the basic variables, and the three equations give their values. At the origin, $x=y=0$ and the basic variables are the slacks: $14,16,6$.

One pivot moves to a neighbour in two choices:

1. Entering variable. Pick a zero variable whose increase raises the profit. At the origin, raising $x$ earns 7 per unit and raising $y$ earns 5. The rule used here, Bland's rule, always takes the smallest-numbered such variable: $x$.
2. Leaving variable. Raise $x$ until some basic variable hits zero. $s_1=14-2x$ hits zero at $x=7$, $s_2=16-x$ at $x=16$, $s_3=6-x$ at $x=6$. The smallest of these ratios, 6, wins: $s_3$ leaves, and we stand at $(6,0)$ with profit 42.

Then rewrite the equations so that $x$ is expressed through the new zero variables (exact row operations on a table of coefficients, the tableau) and repeat. The second pivot raises $y$ and stops where the cutting rule becomes tight, at $(6,2)$ with 52. The third is less obvious: it raises $s_3$, the bench slack, which means making fewer benches. Along the cutting line, trading one bench for two shelves loses 7 and gains 10, so profit rises until the finishing rule becomes tight at $(4,6)$, profit 58. Now every zero variable would lower the profit if raised, and the method stops.

Why the minimum ratio: raise $x$ to 7 instead (leaving by the cutting rule) and $s_3=6-7=-1$; the plan makes 7 benches, and the market takes only 6. Any row other than the minimum-ratio one overshoots some other rule.

Why it stops at the optimum (sketched): when no zero variable can raise the objective, the objective row of the tableau writes the profit as $58$ minus a non-negative combination of the zero variables. Every feasible point has those variables $\ge 0$, so none earns more than 58. Module 16 turns that row into a certificate the partner can check without trusting any code.

Two details make it a real algorithm. When more than $n$ boundaries meet at one corner (a degenerate corner), a pivot can have ratio 0: the basis changes, the point and the value do not. Hence "has not decreased" in the invariant. A run of such pivots can cycle forever under an unlucky rule; Bland's rule (smallest index enters, and among tied ratios the smallest index leaves) provably never repeats a basis, so it terminates. And if the entering variable can grow without any basic variable reaching zero, the LP is unbounded, which is the method's certificate for the second of the three outcomes.

The implementation

This is the course's exact simplex method, which every lab of modules 14–17 provides. pivot does the row operations, run_simplex applies Bland's rule, and simplex builds the tableau with one slack per row. When some $b_i<0$ the origin is not feasible, so simplex first runs a phase 1: it adds one artificial variable $x_0$ and pivots to make $-x_0$ as large as possible; if that maximum is below zero, the LP is infeasible. The on_pivot hook is called after every pivot of the main phase, and it is how we watch the walk.

from fractions import Fraction as F


def pivot(T, B, r, s):
    """Make column s basic in row r (exact row operations on the tableau)."""
    T[r] = [v / T[r][s] for v in T[r]]
    for i in range(len(T)):
        if i != r and T[i][s] != 0:
            f = T[i][s]
            T[i] = [a - f * b for a, b in zip(T[i], T[r])]
    B[r] = s


def run_simplex(T, B, cols, on_pivot=None):
    """Bland's rule: smallest entering index, then smallest leaving index among ties."""
    while True:
        s = next((j for j in cols if T[-1][j] < 0), None)
        if s is None:
            return "optimal"
        rows = [i for i in range(len(B)) if T[i][s] > 0]
        if not rows:
            return "unbounded"
        r = min(rows, key=lambda i: (T[i][-1] / T[i][s], B[i]))
        pivot(T, B, r, s)
        if on_pivot:
            on_pivot(T, B, r, s)


def simplex(c, A, b, on_pivot=None):
    """max c.x s.t. A x <= b, x >= 0, exactly. Returns (status, x, value, y);
    y are the dual prices (y >= 0, A^T y >= c, b.y == value at optimum)."""
    m, n = len(A), len(c)
    art = n + m                                   # auxiliary column x0, used only if some b < 0
    T = [[F(v) for v in A[i]] + [F(int(i == k)) for k in range(m)] + [F(-1), F(b[i])]
         for i in range(m)]
    B = [n + i for i in range(m)]
    if m and min(b) < 0:                          # phase 1: maximize -x0
        T.append([F(0)] * art + [F(1), F(0)])
        pivot(T, B, min(range(m), key=lambda i: b[i]), art)
        run_simplex(T, B, range(art + 1))
        if T[-1][-1] < 0:
            return "infeasible", None, None, None
        if art in B:                              # x0 basic at 0: pivot it out
            r = B.index(art)
            pivot(T, B, r, next(j for j in range(art) if T[r][j] != 0))
        T.pop()
    z = [F(-v) for v in c] + [F(0)] * (m + 2)
    for i, j in enumerate(B):
        if z[j] != 0:
            f = z[j]
            z = [a - f * t for a, t in zip(z, T[i])]
    T.append(z)
    if run_simplex(T, B, range(art), on_pivot) == "unbounded":
        return "unbounded", None, None, None
    x = [F(0)] * (n + m)
    for i, j in enumerate(B):
        x[j] = T[i][-1]
    return "optimal", x[:n], T[-1][-1], [T[-1][n + i] for i in range(m)]

walk = []
def watch(T, B, r, s):
    x = [F(0)] * 5
    for i, j in enumerate(B):
        x[j] = T[i][-1]
    walk.append((s, (x[0], x[1]), T[-1][-1]))

status, plan, profit, _ = simplex([7, 5], [[2, 1], [1, 2], [1, 0]],
                                  [14, 16, 6], watch)
assert (status, plan, profit) == ("optimal", [4, 6], 58)
assert walk == [(0, (6, 0), 42), (1, (6, 2), 52), (4, (4, 6), 58)]
assert all(a[2] <= b[2] for a, b in zip(walk, walk[1:]))   # never downhill

walk = []                     # the other way round: number y before x
simplex([5, 7], [[1, 2], [2, 1], [0, 1]], [14, 16, 6], watch)
assert [(v[::-1], z) for _, v, z in walk] == [((0, 8), 40), ((4, 6), 58)]

Standard form: four rewrites

simplex accepts only $max$, $\le$ and $x\ge 0$. Every LP can be put in that form without changing its answer:

• min to max: minimizing $c·x$ is maximizing $(-c)·x$; negate the optimum back.
• $\ge$ to $\le$: $a·x\ge \beta$ is $(-a)·x\le -\beta$.
• $=$ to two rows: $a·x=\beta$ is $a·x\le \beta$ together with $(-a)·x\le -\beta$.
• a free variable: if $x_j$ may be negative, write $x_j=x_j^{+}-x_j^{-}$ with both parts $\ge 0$, and replace $x_j$ everywhere by the difference.

Lemma. Each rewrite maps feasible points of the old LP to feasible points of the new one with the same objective value, and back. Proof for the free variable: any real $x_j$ is a difference $max(x_j,0)-max(-x_j,0)$ of non-negative numbers, and any pair $(x_j^{+},x_j^{-})$ gives back the real number $x_j^{+}-x_j^{-}$. The constraints and the objective see only the difference, so feasibility and value carry over both ways. That many pairs give the same $x_j$ is harmless. The other three rewrites are one-line identities. $◻$

A diet problem shows the first two. Feed A costs 3 per bag and feed B costs 2. A mix needs at least 9 units of protein ($3a+b\ge 9$) and at least 8 of fibre ($a+2b\ge 8$). Minimize the cost:

status, mix, value, _ = simplex([-3, -2], [[-3, -1], [-1, -2]], [-9, -8])
assert status == "optimal" and mix == [2, 3] and -value == 12

Both rows had negative right-hand sides after the rewrite, so this call ran phase 1. The cheapest mix is 2 bags of A and 3 of B, cost 12.

One more device turns a maximum into linear rules. To minimize the largest of several linear terms, add a variable $z$, require $z\ge$ each term, and minimize $z$. At the optimum, $z$ equals the largest term, since otherwise it could be lowered. Module 14's game value was the mirror image: the row player maximizes $v$ subject to $v\le$ each column's expected payoff.

Modelling: knapsack, flow, shortest path

The hard part of using LPs is choosing variables whose rules are linear. Three patterns cover a great deal.

Fractional knapsack. Items worth 12, 10 and 7 weigh 4, 5 and 7, the bag holds 10, and any fraction of an item may be taken. Let $x_i\in [0,1]$ be the fraction of item $i$: maximize $12x_1+10x_2+7x_3$ subject to $4x_1+5x_2+7x_3\le 10$ and each $x_i\le 1$. The optimum takes the first two whole and $1/7$ of the third, value 23, which is the greedy-by-value-per-weight answer: here the LP confirms a greedy rule.

status, take, value, _ = simplex([12, 10, 7],
                                 [[4, 5, 7], [1, 0, 0], [0, 1, 0], [0, 0, 1]],
                                 [10, 1, 1, 1])
assert take == [1, 1, F(1, 7)] and value == 23

Maximum flow. Give each edge a variable $f_e$ with $0\le f_e\le$ its capacity; at every vertex other than $s$ and $t$, flow in equals flow out (an equality, so two rows); maximize the net flow out of $s$. The feasible points are exactly the flows of module 11, so the optimum is the maximum flow. On the network below the LP's optimum is 9.

NET = {("s", "a"): 6, ("s", "b"): 4, ("a", "b"): 2, ("a", "c"): 3,
       ("b", "c"): 1, ("b", "t"): 5, ("c", "t"): 7}
FLOW = {("s", "a"): 5, ("s", "b"): 4, ("a", "b"): 2, ("a", "c"): 3,
        ("b", "c"): 1, ("b", "t"): 5, ("c", "t"): 4}
assert all(0 <= FLOW[e] <= NET[e] for e in NET)
for v in "abc":
    inflow = sum(f for (p, q), f in FLOW.items() if q == v)
    assert inflow == sum(f for (p, q), f in FLOW.items() if p == v)
assert FLOW[("s", "a")] + FLOW[("s", "b")] == 9
cut = [c for (p, q), c in NET.items() if p in "sa" and q not in "sa"]
assert sum(cut) == 9                               # the cut {s, a}

The flow above is one optimal vertex, and the cut $\{s,a\}$ with capacity $4+2+3=9$ proves, by module 11, that no flow does better. That the LP produces such a certificate on its own is the content of module 16.

Shortest path. On the graph $s\to a$ 4, $s\to b$ 1, $b\to a$ 2, $a\to t$ 5, $b\to t$ 8, give each vertex a free variable $d_v$ with $d_s=0$, and for each edge $u\to v$ of length $w$ require $d_v-d_u\le w$. Then maximize $d_t$.

Lemma. With non-negative lengths and $t$ reachable, the maximum of $d_t$ is the distance from $s$ to $t$. Proof. Along any path $s=v_0,v_1,\dots ,v_k=t$, adding the constraints $d_{v_{i+1}}-d_{v_i}\le w_i$ telescopes to $d_t-d_s\le$ the path's length, so $d_t$ is at most the shortest length. And the true distances satisfy every constraint (the triangle inequality), with $d_t$ equal to the distance. $◻$

Here the distance is 8, along $s\to b\to a\to t$ ($1+2+5$).

EDGES = [("s", "a", 4), ("s", "b", 1), ("b", "a", 2),
         ("a", "t", 5), ("b", "t", 8)]
def fits(d):
    return d["s"] == 0 and all(d[v] - d[u] <= w for u, v, w in EDGES)
assert fits({"s": 0, "a": 3, "b": 1, "t": 8})         # d_t = 8 fits
assert 1 + 2 + 5 == 8                  # the path s b a t caps d_t at 8
assert fits({"s": 0, "a": 3, "b": 1, "t": -10 ** 6})  # no floor

The direction matters. The constraints only bound $d$ from above, so "minimize $d_t$" is unbounded: $d_t=-1,000,000$ breaks no rule. The distance is the largest $d_t$ allowed. The LP also stays correct with negative lengths, as long as no cycle has negative total length (with one, it is infeasible).

How many pivots?

Each pivot rewrites the tableau: $m+1$ rows of $n+m+2$ entries, so $O(m(m+n))$ exact arithmetic operations. The question is how many pivots. The number of bases bounds it, which is as bad as vertex enumeration, but in practice far fewer are needed. Here is a measurement on random LPs with $n$ variables and $n$ rules (entries of $A$ and $c$ in 1..9, $b$ in 10..60, five per size, seed 451):

import random

def make_lp(n, r):
    A = [[r.randint(1, 9) for _ in range(n)] for _ in range(n)]
    b = [r.randint(10, 60) for _ in range(n)]
    c = [r.randint(1, 9) for _ in range(n)]
    return c, A, b

rng, table = random.Random(451), []
for n in (2, 3, 4, 5, 6):
    counts = []
    for _ in range(5):
        pivots = []
        simplex(*make_lp(n, rng), lambda *args: pivots.append(1))
        counts.append(len(pivots))
    table.append((n, comb(2 * n, n), max(counts), sum(counts) / 5))
assert [row[1] for row in table] == [6, 20, 70, 252, 924]
assert [row[2] for row in table] == [2, 2, 4, 6, 15]
assert all(row[2] <= 3 * row[0] for row in table)
assert [round(row[3], 1) for row in table] == [2.0, 1.4, 3.0, 2.8, 6.2]

On these instances simplex never needed more than $3n$ pivots, while enumeration's count grows like $4^n$. That is a measurement on 25 random LPs, not a theorem. No pivot rule is known to need only polynomially many pivots on every LP. For the classic largest-coefficient rule there are LPs that force exponentially many (Klee and Minty, 1972), and Bland's rule has such LPs too; module 17 builds one. LPs can be solved in time polynomial in the number of bits of the input, by the ellipsoid method (Khachiyan, 1979) and by interior-point methods (Karmarkar, 1984). Those bounds are stated here, not proved. Exact fractions also have a price: their numerators and denominators can grow, which is why production solvers use floating point with tolerances.

A problem that looks different

A paint shop mixes three base paints, each with a known share of red, yellow and blue pigment and a price per litre, to make 20 litres whose pigment shares must each land within a tolerance of a customer's sample. It wants the cheapest mix. What would the unknowns be, and is every rule a sum of unknowns times fixed numbers? The lab's last problem is a different one.

Practise

In the lab you list every corner of a new two-variable LP by solving each candidate system, and watch each candidate accepted or rejected; predict the simplex walk on that LP, including a pivot that changes nothing, and then trace it through the hook; write the standard-form front end and use it to compute maximum flows and all shortest-path distances at once, checked against the real algorithms; count pivots on larger random LPs; and solve a planning problem that does not say what it is.

Recap

You can now: draw a two-variable LP, name its corners and find the optimum by sliding the objective line; convert any LP to $maxc·x$, $Ax\le b$, $x\ge 0$; model production, diet, knapsack, maximum flow and shortest path as LPs and say what each variable means; and trace the simplex method pivot by pivot.

Invariant: every simplex pivot moves from a vertex to a neighbouring vertex (or stays, at a degenerate corner) without lowering the objective, and each standard-form rewrite keeps the feasible points and their values.

Complexity achieved: a few pivots of $O(m(m+n))$ operations each on the measured instances (at most $3n$ there), against $\left(\genfrac{}{}{0ex}{}{m+n}{n}\right)$ linear systems for vertex enumeration. There is no known polynomial bound on pivots; the ellipsoid and interior-point methods are polynomial.

Failure mode: a wrong direction (minimizing $d_t$) or a missing sign rule (a free variable left $\ge 0$) silently changes the problem. The solver answers the LP you wrote.

In real software: SciPy's scipy.optimize.linprog uses the HiGHS solvers by default since SciPy 1.9 (method='highs', choosing between a dual simplex method and an interior-point method), and by default gives every variable the bounds $(0,\infty )$.

Retrieval (module 11): on this lesson's network, how do you know the flow of 9 is maximum without trusting the LP?

Check yourself

1. Why is the workshop's optimum at a corner, and why at $(4,6)$ rather than $(6,2)$?
2. Raise the shelf profit from 5 to 15. What changes, and at what shelf profit does the optimum first leave $(4,6)$?
3. Solve this lesson's shortest-path instance by Dijkstra's algorithm and by the LP. Which one still works if an edge has negative length, and which is faster?

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.