iv. Linear Programming

Optimise a linear objective subject to linear inequalities:
$\max c^{\top }x\text{s.t.}Ax\le b,x\ge 0.$
General LPs: the simplex method (exponential worst case, a few times $m$ pivots in practice). Fixed dimension $d\le 3$ with many constraints: Seidel’s randomized algorithm in expected $O(d!n)$, or half-plane intersection when the whole region is needed.

Have	Convert to standard form
minimise $c^{\top }x$	maximise $-c^{\top }x$
$a^{\top }x\ge b$	$-a^{\top }x\le -b$
$a^{\top }x=b$	two inequalities $\le b$ and $\ge b$
free variable $x$	$x=x^{+}-x^{-}$ with $x^{\pm }\ge 0$
$|a^{\top }x|\le b$	two inequalities
minimise ${\max }_i{a}_i^{\top }x$	new variable $t$, constraints $a_i^{\top }x\le t$, minimise $t$

Duality

Primal: $\max c^{\top }x$, $Ax\le b$, $x\ge 0$. Dual: $\min b^{\top }y$, $A^{\top }y\ge c$, $y\ge 0$.

• Weak duality: every primal value $\le$ every dual value.
• Strong duality: if the primal has a finite optimum, so does the dual, and they are equal.
• Complementary slackness: at optima, $y_i>0\Rightarrow$ primal constraint $i$ is tight, and $x_j>0\Rightarrow$ dual constraint $j$ is tight.
• Total unimodularity: if every square submatrix of $A$ has determinant in $\{-1,0,1\}$ (network matrices, bipartite incidence matrices, interval matrices with consecutive ones) and $b$ is integral, every vertex of $\{Ax\le b,x\ge 0\}$ is integral: the LP relaxation solves the integer problem, usually via flows.

Duality explains max-flow = min-cut and Kőnig’s theorem, gives potentials in min-cost flow, and lets you swap a hard LP for an easier dual: many “minimum cost to satisfy covering constraints” problems are duals of flows (Flow Modeling).

Proof: weak duality, complementary slackness and integrality under TU

Weak duality. For feasible $x,y$: $c^{\top }x\le (A^{\top }y{)}^{\top }x=y^{\top }Ax\le y^{\top }b$, using $x\ge 0$ with $A^{\top }y\ge c$, then $y\ge 0$ with $Ax\le b$. Strong duality is the LP duality theorem (Farkas’ lemma); constructively, when the simplex method stops with all reduced costs non-negative, the multipliers of the final basis form a dual-feasible $y$ with $b^{\top }y=c^{\top }x$. Complementary slackness. At optima equality holds in both inequalities above: ${\sum }_j{x}_j((A^{\top }y{)}_j-c_j)=0$ and ${\sum }_i{y}_i(b_i-(Ax{)}_i)=0$, sums of non-negative terms, so every term is 0. Total unimodularity. A vertex is the unique solution of a square subsystem $A_B{x}_B=b_B$ (tight constraints, including some $x_j=0$). If $A$ is TU, so is $[A;I]$, the subsystem’s determinant is $\pm 1$, and by Cramer’s rule $x_B=A_B^{-1}{b}_B$ is integral.

Simplex (two-phase, Bland’s rule)

LPSolver(A, b, c).solve(x) returns the optimum and fills $x$; it returns $-\infty$ if infeasible and $+\infty$ if unbounded. $O(nm)$ per pivot; fine for $n,m$ up to a few hundred. Phase 1 finds a feasible basis when some $b_i<0$ (the origin is infeasible).

using T = long double; using vd = vector<T>; using vvd = vector<vd>;
const T LP_EPS = 1e-9, LP_INF = numeric_limits<T>::infinity();
 
struct LPSolver {                                          // maximize c x, A x <= b, x >= 0 (KACTL-style)
    int m, n; vector<int> N, B; vvd D;
    LPSolver(const vvd& A, const vd& b, const vd& c) : m(sz(b)), n(sz(c)), N(n + 1), B(m), D(m + 2, vd(n + 2)) {
        for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) D[i][j] = A[i][j];
        for (int i = 0; i < m; i++) { B[i] = n + i; D[i][n] = -1; D[i][n + 1] = b[i]; }
        for (int j = 0; j < n; j++) { N[j] = j; D[m][j] = -c[j]; }
        N[n] = -1; D[m + 1][n] = 1;
    }
    void pivot(int r, int s) {
        T inv = 1 / D[r][s];
        for (int i = 0; i < m + 2; i++) if (i != r && fabsl(D[i][s]) > LP_EPS) {
            T f = D[i][s] * inv;
            for (int j = 0; j < n + 2; j++) D[i][j] -= D[r][j] * f;
            D[i][s] = D[r][s] * f;
        }
        for (int j = 0; j < n + 2; j++) if (j != s) D[r][j] *= inv;
        for (int i = 0; i < m + 2; i++) if (i != r) D[i][s] *= -inv;
        D[r][s] = inv; swap(B[r], N[s]);
    }
    bool simplex(int phase) {
        int x = m + phase - 1;
        while (true) {
            int s = -1;
            for (int j = 0; j <= n; j++) if (N[j] != -phase) if (s == -1 || make_pair(D[x][j], N[j]) < make_pair(D[x][s], N[s])) s = j;
            if (D[x][s] >= -LP_EPS) return true;          // optimal
            int r = -1;
            for (int i = 0; i < m; i++) {
                if (D[i][s] <= LP_EPS) continue;
                if (r == -1 || make_pair(D[i][n + 1] / D[i][s], B[i]) < make_pair(D[r][n + 1] / D[r][s], B[r])) r = i;
            }
            if (r == -1) return false;                     // unbounded
            pivot(r, s);
        }
    }
    T solve(vd& x) {
        int r = 0;
        for (int i = 1; i < m; i++) if (D[i][n + 1] < D[r][n + 1]) r = i;
        if (D[r][n + 1] < -LP_EPS) {                       // phase 1: find a feasible basis
            pivot(r, n);
            if (!simplex(2) || D[m + 1][n + 1] < -LP_EPS) return -LP_INF;   // infeasible
            for (int i = 0; i < m; i++) if (B[i] == -1) {
                int s = 0;
                for (int j = 1; j <= n; j++) if (make_pair(D[i][j], N[j]) < make_pair(D[i][s], N[s])) s = j;
                pivot(i, s);
            }
        }
        bool ok = simplex(1); x = vd(n);
        for (int i = 0; i < m; i++) if (B[i] < n) x[B[i]] = D[i][n + 1];
        return ok ? D[m][n + 1] : LP_INF;
    }
};
// maximize 3x + 2y s.t. x + y <= 4, x + 3y <= 6, x, y >= 0  ->  12 at x = (4, 0)

Pitfalls. Variables are implicitly $\ge 0$ (split free variables); equalities become two inequalities; use long double and $\epsilon =10^{-9}$, and scale constraints to similar magnitudes; degenerate pivots can cycle, which Bland’s rule (smallest index among ties) prevents; distinguish infeasible from unbounded in the output; the answer is floating point: round only if integrality is guaranteed (TU). Integer LP is NP-hard.

Proof: what the simplex code guarantees

The tableau always represents the LP in the current basis: row $i$ expresses basic variable $B_i$ through the non-basic ones, D[m] holds the reduced costs of the objective and D[m][n+1] its value at the basic solution (non-basic variables 0). pivot is Gauss–Jordan elimination exchanging $B_r$ and $N_s$, so the represented LP is unchanged. Optimality: if no reduced cost is negative, increasing any non-basic variable cannot raise the objective, and the current basic solution (feasible: right-hand sides $\ge 0$ are preserved by the ratio test) is optimal; the final multipliers are a dual certificate. Unboundedness: if the entering column has no positive entry, increasing that variable keeps all basic variables non-negative forever while the objective grows. Termination: with Bland’s rule (smallest index among candidates, implemented by the make_pair tie-breaks) the simplex method never cycles (Bland, 1977), and there are finitely many bases. Phase 1: the auxiliary variable (index $-1$, column $n$) is added to every row; pivoting it in on the most negative $b_r$ makes all right-hand sides non-negative, and maximizing $-x_{\text{aux}}$ reaches 0 iff the original system is feasible, after which the auxiliary variable is pivoted out. Checked against brute-force vertex enumeration on random LPs with up to 3 variables and 6 constraints, including infeasible and unbounded instances.

Seidel’s LP in fixed dimension

Shuffle the constraints and add them one at a time, maintaining the current optimum $x^{\ast }$. If $x^{\ast }$ satisfies the new constraint, nothing changes (the common case). Otherwise the new optimum lies on the boundary of the new constraint: project the earlier constraints onto that hyperplane and solve a $(d-1)$-dimensional LP. Expected time $O(d!n)$, space $O(n)$; the deterministic alternative (Megiddo) is $O(2^{2^d}n)$. Below is the 2D version, with a bounding box $|x|,|y|\le M$ that makes every LP bounded.

struct HP { long double a, b, c; };                        // a x + b y <= c
 
// maximize cx * x + cy * y subject to the half-planes and |x|, |y| <= M. Returns {feasible, (x, y)}.
pair<bool, pair<long double, long double>> seidel2D(vector<HP> H, long double cx, long double cy, long double M, mt19937_64& rng) {
    const long double eps = 1e-9;
    shuffle(all(H), rng);
    vector<HP> box = {{1, 0, M}, {-1, 0, M}, {0, 1, M}, {0, -1, M}};
    H.insert(H.begin(), all(box));
    long double x = cx >= 0 ? M : -M, y = cy >= 0 ? M : -M;   // optimum of the box alone
    for (int i = 4; i < sz(H); i++) {
        auto [a, b, c] = H[i];
        if (a * x + b * y <= c + eps) continue;
        long double l2 = a * a + b * b;
        if (l2 < eps) return {false, {0, 0}};               // 0 <= c violated
        long double px = a * c / l2, py = b * c / l2, dx = -b, dy = a;   // the line: (px, py) + t (dx, dy)
        long double lo = -1e30, hi = 1e30;
        for (int j = 0; j < i; j++) {
            long double num = H[j].c - (H[j].a * px + H[j].b * py), den = H[j].a * dx + H[j].b * dy;
            if (fabsl(den) < eps) { if (num < -eps) return {false, {0, 0}}; continue; }
            if (den > 0) hi = min(hi, num / den); else lo = max(lo, num / den);
        }
        if (lo > hi + eps) return {false, {0, 0}};
        long double dir = cx * dx + cy * dy, t = dir > 0 ? hi : lo;
        x = px + t * dx; y = py + t * dy;
    }
    return {true, {x, y}};
}

Backwards analysis is the real lesson: after $i$ constraints the optimum is determined by $d$ of them, and in random order the $i$-th constraint is one of those with probability $\frac{d}{i}$, so $T_d(n)=O(n)+{\sum }_i\frac{d}{i}{T}_{d-1}(i)=O(d!n)$. The same argument gives Welzl’s circle ($\mathrm{Pr}=\frac{3}{i}$), randomized incremental Delaunay ($O(1)$ expected change per insertion), randomized quicksort ($O(n\log n)$) and treap depth ($O(\log n)$). Seidel’s algorithm extends to every LP-type problem (combinatorial dimension $d$, monotonicity and locality): minimum enclosing ball and ellipsoid, smallest enclosing annulus, distance between polytopes (Sharir–Welzl; Matoušek–Sharir–Welzl).

Contest use of low-dimensional LP
is the intersection of $n$ half-planes non-empty?	2D LP feasibility
Chebyshev centre (largest inscribed circle of a convex polygon)	3D LP in $(x,y,r)$
smallest enclosing circle	LP-type: Welzl
line fitting under $L_{\infty }$ (minimise the maximum vertical error)	2D LP
difference constraints	better with Bellman–Ford
”minimise a max of linear functions” in 1 variable	ternary search (convex); 2–3 variables: nested ternary search or LP
fractional relaxations	bounds in branch and bound

Proof: Seidel's algorithm is correct

Invariant: after processing constraint $i$, $(x,y)$ is an optimal point of the LP with the box and constraints $\mathrm{1..}i$. If the current point satisfies constraint $i$, it remains optimal (the feasible region only shrank). Otherwise let $q$ be any optimal point of the new LP (if it is feasible) and $p$ the old optimum, which violates constraint $i$ and has $c\cdot p\ge c\cdot q$. On the segment from $q$ to $p$, the point $z$ where it crosses the boundary line of constraint $i$ satisfies all earlier constraints (convexity) and constraint $i$, and $c\cdot z\ge c\cdot q$ (linear interpolation): so an optimum lies on the line. The 1D LP on the line intersects the earlier constraints into an interval $[lo,hi]$ (finite thanks to the box) and takes the endpoint in the objective’s direction; an empty interval proves infeasibility of the whole LP. Expected time follows from backwards analysis: constraint $i$ forces the $O(i)$ recomputation only if it is one of the (at most 2) constraints defining the optimum after $i$ insertions, probability $\le \frac{2}{i}$, total $O(n)$. Checked against brute-force enumeration of all pairwise line intersections on random instances.

Related notes

• Flow Modeling (LPs that are flows)
• Half-Plane Intersection
• Approximation Algorithms (LP rounding)