Simplex Algorithm, in depth: vertices, pivots, duals, degeneracy and what real LP solvers do

Linear programs and standard form

A linear program in the form we will use is: maximise c·x subject to A x ≤ b and x ≥ 0, where A is an m-by-n matrix. Other shapes convert mechanically. Minimising c·x is maximising -c·x. An equality is two inequalities. A ≥ constraint is a ≤ constraint with both sides negated. A free variable is the difference of two non-negative ones.

The method works on equalities, so each inequality gets a slack variable: x1 ≤ 4 becomes x1 + s1 = 4 with s1 ≥ 0. The slack measures unused capacity. Now there are m equations in n + m variables. Pick any m variables whose columns are linearly independent, set the other n to zero, and solve: that is a basic solution. The m chosen variables are the basis; the rest are non-basic. If every value comes out non-negative, it is a basic feasible solution.

The running example is the textbook Wyndor Glass problem, two products sharing three plants' hours: maximise 3x1 + 5x2 subject to x1 ≤ 4, 2x2 ≤ 12, 3x1 + 2x2 ≤ 18, x1, x2 ≥ 0.

The geometry: why the answer is a corner

Each inequality cuts space in half, so the feasible region is an intersection of half-spaces: a convex polyhedron. A linear objective has no curvature, so if a maximum exists on a bounded region, one is attained at a vertex. Sweep the level line 3x1 + 5x2 = z outward and the last point of the polygon it touches is a corner (or an edge whose two corners tie).

Two facts make this an algorithm. First, vertices are exactly the basic feasible solutions. Second, convexity makes local optimality global. If no edge leaving a vertex improves the objective, no point of the polyhedron does. So the simplex method walks from vertex to adjacent vertex along improving edges, and stops when there is none.

Worked example: two pivots by hand

Write the problem as a tableau. Each row is one equation; the bottom row holds the objective written as z - 3x1 - 5x2 = 0. Start at the origin, which is feasible because b is non-negative: the slacks are basic, s1 = 4, s2 = 12, s3 = 18, and z = 0.

basis |  x1   x2   s1   s2   s3 |  rhs
 s1   |   1    0    1    0    0 |    4
 s2   |   0    2    0    1    0 |   12
 s3   |   3    2    0    0    1 |   18
  z   |  -3   -5    0    0    0 |    0

Choose the entering variable. A negative entry in the z row means raising that variable raises z. Dantzig's original rule takes the most negative: x2, worth 5 per unit.

Choose the leaving variable with the ratio test. As x2 grows, row s2 hits zero at 12 / 2 = 6 and row s3 at 18 / 2 = 9. The smallest ratio decides: s2 leaves, because a larger step would push it negative.

Pivot. Divide the s2 row by 2 so x2 has coefficient 1, then subtract multiples of it from every other row, including z, so x2's column becomes a unit vector.

basis |  x1   x2   s1   s2   s3 |  rhs
 s1   |   1    0    1    0    0 |    4
 x2   |   0    1    0   1/2   0 |    6
 s3   |   3    0    0   -1    1 |    6
  z   |  -3    0    0   5/2   0 |   30

We are at vertex (0, 6) with z = 30. x1 still has a negative entry, so it enters. Ratios: row s1 gives 4 / 1 = 4, row s3 gives 6 / 3 = 2. s3 leaves. After the second pivot the z row reads 0 0 0 3/2 1 | 36. No entry is negative, so the vertex x1 = 2, x2 = 6 is optimal with z = 36. Check: 3·2 + 5·6 = 36, and plants 2 and 3 are fully used while plant 1 has 2 hours spare.

A complete implementation

The implementation follows the tableau exactly, using Fraction so it pivots without rounding and serves as a reference for testing a floating-point version.

from fractions import Fraction as F


def simplex_max(c, A, b, rule="bland", max_iter=10_000):
    """Maximise c.x subject to A x <= b, x >= 0, b >= 0. Exact arithmetic."""
    m, n = len(A), len(c)
    # Tableau rows: [A | I | b]; objective row holds reduced costs -c.
    T = [[F(v) for v in A[i]] + [F(int(i == j)) for j in range(m)] + [F(b[i])]
         for i in range(m)]
    obj = [F(-v) for v in c] + [F(0)] * m + [F(0)]
    basis = [n + i for i in range(m)]          # slacks start basic

    for _ in range(max_iter):
        candidates = [j for j in range(n + m) if obj[j] < 0]
        if not candidates:
            break                              # optimal: no improving column
        if rule == "bland":
            e = min(candidates)                # smallest index: cannot cycle
        else:
            e = min(candidates, key=lambda j: obj[j])   # Dantzig: most negative
        ratios = [(T[i][-1] / T[i][e], basis[i], i) for i in range(m) if T[i][e] > 0]
        if not ratios:
            return "unbounded", None, None, None
        _, _, r = min(ratios)                  # ties broken by smallest basic index
        piv = T[r][e]
        T[r] = [v / piv for v in T[r]]
        for i in range(m):
            if i != r and T[i][e] != 0:
                f = T[i][e]
                T[i] = [a - f * p for a, p in zip(T[i], T[r])]
        f = obj[e]
        obj = [a - f * p for a, p in zip(obj, T[r])]
        basis[r] = e
    else:
        return "iteration_limit", None, None, None

    x = [F(0)] * (n + m)
    for i, j in enumerate(basis):
        x[j] = T[i][-1]
    duals = obj[n:n + m]                       # reduced costs of the slacks
    return "optimal", x[:n], obj[-1], duals


print(simplex_max([3, 5], [[1, 0], [0, 2], [3, 2]], [4, 12, 18]))
# optimal, x = [2, 6], z = 36, duals = [0, 3/2, 1]
print(simplex_max([1, 1], [[1, -1]], [1]))     # ('unbounded', None, None, None)

Unboundedness shows up in the ratio test: if the entering column has no positive entry, increasing the variable never makes any basic variable hit zero, so z grows without limit. In the second call, x2 can rise forever. In a real model that almost always means a missing constraint.

Reading the final tableau: duals and reduced costs

The entries under the slack columns, 0, 3/2 and 1, are the dual values or shadow prices: the rate at which the optimum rises per extra unit of each resource. One more hour in plant 2 is worth 1.5; one more in plant 3 is worth 1; plant 1 has spare hours, so another is worth nothing. Strong duality says these prices, multiplied by the capacities, reproduce the optimum: 0·4 + 1.5·12 + 1·18 = 36.

Shadow prices hold only while the basis stays optimal, a range solvers report. The entries under the original variables are reduced costs: how much a non-basic variable's profit would have to rise before it is worth producing. The dual itself is a linear program, minimising b·y subject to Aᵀy ≥ c, y ≥ 0, and the max-flow min-cut theorem is one famous instance of this duality.

Finding a starting vertex: the two-phase method

Starting at the origin only works when b ≥ 0 and every constraint is ≤. A constraint such as x1 + x2 ≥ 2 makes the origin infeasible, and there is no obvious basis. The standard fix is the two-phase method. Phase 1 adds an artificial variable to each awkward row (x1 + x2 - s + a = 2) so the artificials form an initial basis, then minimises their sum with the same pivoting. A positive minimum proves the problem infeasible; a zero minimum yields a feasible vertex from which phase 2 optimises the real objective.

Degeneracy, cycling and pivot rules

A vertex is degenerate when more constraints than necessary are tight there, which appears as a basic variable equal to zero. The ratio test then returns 0, so the pivot changes the basis without moving: z stays put. Real models, especially network-shaped ones, are heavily degenerate.

Zero-length pivots can return to a visited basis and loop forever; this cycling is rare but real. Robert Bland showed in 1977 that choosing the entering variable with the smallest index among improving ones, and breaking ratio-test ties by smallest index too, can never cycle. The code above defaults to it; the price is speed, three pivots instead of two on the running example. Production solvers prevent stalling differently, with small perturbations of b or shifted bounds that are removed at the end.

How many pivots?

Each pivot costs O(m(n + m)) on a dense tableau. Klee and Minty (1972) built squashed cubes on which Dantzig's rule visits all 2ⁿ vertices, so the worst case is exponential, yet in practice the count is usually a small multiple of m. Spielman and Teng's smoothed analysis explains the gap: under small random perturbations of the data, the expected pivot count of a shadow-vertex rule is polynomial.

Polynomial-time LP algorithms exist, the ellipsoid method (Khachiyan, 1979) and interior-point methods (Karmarkar, 1984), which cross the inside of the polyhedron in a few expensive iterations.

What production solvers do

Nobody runs a dense tableau at scale. Production codes use the revised simplex method: keep only the basis matrix B in factored form (a sparse LU factorisation updated after each pivot and refactored periodically), compute the duals y = c_B B⁻¹, price the non-basic columns to pick an entering variable, and compute only the column needed for the ratio test. Pricing uses steepest-edge or Devex weights rather than the most negative entry.

The dual simplex method runs the same machinery on the dual: it keeps the objective row optimal and works towards feasibility. It is the workhorse inside branch and bound, because adding a cut or tightening a bound leaves the previous basis dual feasible, so the solver re-optimises in a few pivots instead of starting over. See Integer Programming and LP, in depth for how that warm start drives an integer solver. Presolve and scaling surround all of this.

In Python, SciPy's linprog wraps the HiGHS solver. The legacy 'simplex', 'revised simplex' and 'interior-point' methods are deprecated (SciPy 1.17 still runs them with a warning); use 'highs-ds' (dual revised simplex), 'highs-ipm' (interior point) or 'highs' to let it choose.

from scipy.optimize import linprog

# linprog minimises, so negate the profits.
res = linprog(c=[-3, -5],
              A_ub=[[1, 0], [0, 2], [3, 2]],
              b_ub=[4, 12, 18],
              bounds=[(0, None), (0, None)],
              method="highs-ds")
print(res.status, res.x, -res.fun)        # 0 [2. 6.] 36.0
print(res.ineqlin.marginals)              # [-0.  -1.5 -1. ]  shadow prices, sign flipped

The marginals are negated because linprog minimises. Check res.status before reading res.x: 2 means infeasible, 3 unbounded.

Failure modes

• Unbounded by accident. A missing upper bound or a sign error produces an unbounded status or an absurd optimum. Put explicit bounds on every variable that has a physical limit.
• Infeasible with no explanation. Phase 1 proves infeasibility but not which constraints conflict. Ask the solver for an irreducible infeasible subsystem where supported, or add penalised slacks.
• Numerical trouble. Coefficients spanning many orders of magnitude, such as 1e-6 and 1e6 in one row, make the basis ill-conditioned and the answer slightly infeasible. Rescale units.
• Float comparisons without tolerance. In a float version, obj[j] < 0 must become obj[j] < -1e-9 or round-off triggers useless pivots.
• Trusting duals outside their range. A shadow price is a local derivative, not a price for 100 extra hours.

Trade-offs

Flow-shaped models have faster dedicated algorithms: see Advanced Flow Networks, in depth and bipartite matching.

What to do next

• Run simplex_max on it, then compare against linprog(method='highs-ds') on random small LPs with non-negative b.
• Add a phase 1 to the implementation and test it on a problem with a ≥ constraint and on one that is infeasible.
• Read the duals of a model you care about and write down, in business terms, what each shadow price means.
• Model a real planning problem with explicit bounds on every variable, solve it with HiGHS, and inspect the status, the duals and the reduced costs.