2 Linear Programming: Vertices, Not Slopes

Chapter 1’s whole toolkit — gradients, curvature, line search — assumes f is smooth. A linear program throws that out: the objective is linear (no curvature at all, so Newton’s step is undefined) and the feasible region is carved out by linear constraints, not free $\mathbb{R}^n$. The good news is that this makes the problem easier in a specific, geometric way, not harder.

2.1 Standard form

Definition — Linear program

\[ \min_{x} \; c^\top x \quad \text{s.t.} \quad A_{eq} x = b_{eq}, \quad A_{ub} x \le b_{ub}, \quad x \ge 0. \]

optimlab.optimizers.linear_programming.LinearProgram holds exactly these five pieces (c, A_eq, b_eq, A_ub, b_ub); A_eq/A_ub may each be omitted if a problem doesn’t need them.

2.2 Why the optimum sits at a vertex

Theorem — LP optima are vertices

If a linear program has an optimal solution, at least one optimal solution is a vertex (extreme point) of the feasible polytope.

A vertex is a point where enough constraints are tight (active) that no feasible direction stays feasible in both signs — informally, “a corner.” The proof is short: walk from any feasible point along a direction that doesn’t change the objective (one exists unless you’re already at a vertex) until you hit a new constraint boundary; repeat until stuck, which must happen at a vertex. This is why simplex — a search that only ever visits vertices — is a complete algorithm rather than a heuristic: it isn’t skipping the interior because interior points are hard to evaluate, it’s skipping them because they’re provably never where the answer is.

2.3 Watching it happen: the feasible region and the path across it

Code

from optimlab.optimizers.linear_programming import LinearProgram, simplex
from optimlab.viz import polytope_figure

# maximize 3x0 + 5x1 s.t. x0<=4, 2x1<=12, 3x0+2x1<=18, x>=0 -- classic textbook LP,
# written here as minimize -3x0-5x1 (optimlab.optimizers.linear_programming always minimizes).
lp = LinearProgram(name="classic_lp", c=[-3.0, -5.0], A_ub=[[1, 0], [0, 2], [3, 2]], b_ub=[4, 12, 18])
result = simplex(lp)
polytope_figure(lp, result).show()

Figure 2.1: The feasible region (shaded) and the exact sequence of vertices simplex visited to reach the optimum (star).

Simplex starts at the origin (always feasible here, since every constraint is <= with a nonnegative right-hand side) and moves to an adjacent vertex only if doing so improves the objective — result.n_iter pivots, result.vertices the exact path drawn above. tests/test_linear_programming.py::test_records_a_vertex_per_pivot checks every recorded point really is a vertex of the true feasible region, not an intermediate approximation.

2.4 How simplex decides which vertex is next

Each pivot picks an entering variable — a nonbasic variable whose reduced cost is negative, meaning increasing it from 0 would improve the objective — and a leaving variable, chosen by a ratio test that finds how far the entering variable can increase before some currently-basic variable would go negative. optimlab.optimizers .linear_programming picks both by Bland’s rule (smallest eligible index, not Dantzig’s largest-improvement rule) specifically because it has a clean proof of termination: with an arbitrary tie-breaking rule, simplex can in principle cycle forever between a set of degenerate vertices; Bland’s rule provably can’t, at the cost of sometimes taking a less direct path than “always pick the most negative reduced cost” would.

Most linear programs don’t start with an obvious feasible vertex to walk from — a problem built entirely from equality constraints has no free slack variables at all. simplex handles this with the standard two-phase construction: phase 1 minimizes the sum of temporary artificial variables (a placeholder that’s easy to start feasible, e.g. one per equality row) starting from an easy-to-write-down infeasible point; reaching an artificial-sum of exactly zero is a feasible vertex of the real problem, which phase 2 then optimizes from. tests/test_linear_programming.py ::test_equality_constraints_need_phase_one exercises exactly this path.

2.5 Correctness: two independent solvers agree

Code

from optimlab.backends import scipy_linprog

oracle = scipy_linprog(lp)
print(f"ours:  x={result.x}, objective={result.objective:.6g}")
print(f"scipy: x={oracle.x}, objective={oracle.objective:.6g}")

ours:  x=[2. 6.], objective=-36
scipy: x=[2. 6.], objective=-36

optimlab.backends.scipy_linprog wraps scipy’s HiGHS solver — an entirely different implementation (dual simplex plus presolve, not our from-scratch primal tableau) landing on the same vertex is meaningfully stronger evidence than either solver’s self-report. optimlab.backends.cvxpy_linprog (needs the backends extra) provides a third, independent check via a different solver stack again (Clarabel via cvxpy); see tests/test_backends_cvxpy.py for a case where scipy and cvxpy were cross-checked against each other directly.

2.6 Infeasible and unbounded problems

Not every LP has an optimum. simplex reports which failure mode explicitly rather than returning a nonsense point:

Code

infeasible = LinearProgram(c=[1.0, 1.0], A_ub=[[1, 1], [-1, -1]], b_ub=[1, -3])  # x0+x1<=1 AND >=3
unbounded = LinearProgram(c=[-1.0, -1.0], A_ub=[[1.0, -1.0]], b_ub=[1.0])  # x1 can grow forever

print("infeasible problem status:", simplex(infeasible).status)
print("unbounded problem status:", simplex(unbounded).status)

infeasible problem status: infeasible
unbounded problem status: unbounded

Phase 1 is what catches infeasibility: if the minimum possible sum of artificial variables is greater than zero, no point exists that satisfies every original constraint — see the status != "optimal" or obj_row[-1] < -1e-6 check in simplex’s source. Unboundedness is caught during an ordinary pivot: if the entering variable’s column has no positive entry anywhere, increasing it forever never violates any constraint, so the objective can improve forever too.

2.7 What’s next

Chapter 3 stays inside convexity but changes the objective from linear to quadratic (least squares) — where, unlike here, the optimum usually sits in the interior, not at a vertex, and the relevant question shifts from “which vertex” to “how sensitive is the answer to noise in the data,” via the singular value decomposition.

--- title: "Linear Programming: Vertices, Not Slopes" --- Chapter 1's whole toolkit — gradients, curvature, line search — assumes `f` is smooth. A linear program throws that out: the objective is linear (no curvature at all, so Newton's step is undefined) and the feasible region is carved out by linear constraints, not free $\mathbb{R}^n$. The good news is that this makes the problem *easier* in a specific, geometric way, not harder. ## Standard form ::: {.callout-note title="Definition — Linear program"} $$ \min_{x} \; c^\top x \quad \text{s.t.} \quad A_{eq} x = b_{eq}, \quad A_{ub} x \le b_{ub}, \quad x \ge 0. $$ ::: `optimlab.optimizers.linear_programming.LinearProgram` holds exactly these five pieces (`c`, `A_eq`, `b_eq`, `A_ub`, `b_ub`); `A_eq`/`A_ub` may each be omitted if a problem doesn't need them. ## Why the optimum sits at a vertex ::: {.callout-tip title="Theorem — LP optima are vertices"} If a linear program has an optimal solution, at least one optimal solution is a vertex (extreme point) of the feasible polytope. ::: A vertex is a point where enough constraints are tight (active) that no feasible direction stays feasible in both signs — informally, "a corner." The proof is short: walk from any feasible point along a direction that doesn't change the objective (one exists unless you're already at a vertex) until you hit a new constraint boundary; repeat until stuck, which must happen at a vertex. This is *why* simplex — a search that only ever visits vertices — is a complete algorithm rather than a heuristic: it isn't skipping the interior because interior points are hard to evaluate, it's skipping them because they're provably never where the answer is. ## Watching it happen: the feasible region and the path across it ```{python} #| label: fig-lp-polytope #| fig-cap: "The feasible region (shaded) and the exact sequence of vertices simplex visited to reach the optimum (star)." from optimlab.optimizers.linear_programming import LinearProgram, simplex from optimlab.viz import polytope_figure # maximize 3x0 + 5x1 s.t. x0<=4, 2x1<=12, 3x0+2x1<=18, x>=0 -- classic textbook LP, # written here as minimize -3x0-5x1 (optimlab.optimizers.linear_programming always minimizes). lp = LinearProgram(name="classic_lp", c=[-3.0, -5.0], A_ub=[[1, 0], [0, 2], [3, 2]], b_ub=[4, 12, 18]) result = simplex(lp) polytope_figure(lp, result).show() ``` Simplex starts at the origin (always feasible here, since every constraint is `<=` with a nonnegative right-hand side) and moves to an *adjacent* vertex only if doing so improves the objective — `result.n_iter` pivots, `result.vertices` the exact path drawn above. `tests/test_linear_programming.py::test_records_a_vertex_per_pivot` checks every recorded point really is a vertex of the true feasible region, not an intermediate approximation. ## How simplex decides which vertex is next Each pivot picks an **entering variable** — a nonbasic variable whose reduced cost is negative, meaning increasing it from 0 would improve the objective — and a **leaving variable**, chosen by a ratio test that finds how far the entering variable can increase before some currently-basic variable would go negative. `optimlab.optimizers .linear_programming` picks both by **Bland's rule** (smallest eligible index, not Dantzig's largest-improvement rule) specifically because it has a clean proof of termination: with an arbitrary tie-breaking rule, simplex can in principle cycle forever between a set of degenerate vertices; Bland's rule provably can't, at the cost of sometimes taking a less direct path than "always pick the most negative reduced cost" would. Most linear programs don't start with an obvious feasible vertex to walk from — a problem built entirely from equality constraints has no free slack variables at all. `simplex` handles this with the standard **two-phase** construction: phase 1 minimizes the sum of temporary *artificial* variables (a placeholder that's easy to start feasible, e.g. one per equality row) starting from an easy-to-write-down infeasible point; reaching an artificial-sum of exactly zero *is* a feasible vertex of the real problem, which phase 2 then optimizes from. `tests/test_linear_programming.py ::test_equality_constraints_need_phase_one` exercises exactly this path. ## Correctness: two independent solvers agree ```{python} #| label: lp-oracle-check from optimlab.backends import scipy_linprog oracle = scipy_linprog(lp) print(f"ours: x={result.x}, objective={result.objective:.6g}") print(f"scipy: x={oracle.x}, objective={oracle.objective:.6g}") ``` `optimlab.backends.scipy_linprog` wraps scipy's HiGHS solver — an entirely different implementation (dual simplex plus presolve, not our from-scratch primal tableau) landing on the same vertex is meaningfully stronger evidence than either solver's self-report. `optimlab.backends.cvxpy_linprog` (needs the `backends` extra) provides a third, independent check via a different solver stack again (Clarabel via cvxpy); see `tests/test_backends_cvxpy.py` for a case where scipy and cvxpy were cross-checked against each other directly. ## Infeasible and unbounded problems Not every LP has an optimum. `simplex` reports which failure mode explicitly rather than returning a nonsense point: ```{python} #| label: lp-infeasible-unbounded infeasible = LinearProgram(c=[1.0, 1.0], A_ub=[[1, 1], [-1, -1]], b_ub=[1, -3]) # x0+x1<=1 AND >=3 unbounded = LinearProgram(c=[-1.0, -1.0], A_ub=[[1.0, -1.0]], b_ub=[1.0]) # x1 can grow forever print("infeasible problem status:", simplex(infeasible).status) print("unbounded problem status:", simplex(unbounded).status) ``` Phase 1 is what catches infeasibility: if the *minimum possible* sum of artificial variables is greater than zero, no point exists that satisfies every original constraint — see the `status != "optimal" or obj_row[-1] < -1e-6` check in `simplex`'s source. Unboundedness is caught during an ordinary pivot: if the entering variable's column has no positive entry anywhere, increasing it forever never violates any constraint, so the objective can improve forever too. ## What's next Chapter 3 stays inside convexity but changes the objective from linear to quadratic (least squares) — where, unlike here, the optimum usually sits in the *interior*, not at a vertex, and the relevant question shifts from "which vertex" to "how sensitive is the answer to noise in the data," via the singular value decomposition.