I want to explain linear programming (LP) as a mathematical framework I use to find the best outcome — maximum profit, minimum cost, and so on — from a linear objective function, subject to a set of linear constraints. I consider LP one of the most practically important tools in optimization, since a huge range of real-world resource-allocation problems can be modeled as linear programs and solved efficiently, even when they involve thousands or millions of variables.
History and Background
I credit the formal development of linear programming largely to two independent lines of work. Leonid Kantorovich developed early linear-programming-style methods in the Soviet Union in 1939 for production planning. Independently, George Dantzig developed the simplex method in 1947 while working for the U.S. Air Force on logistics and planning problems, and his work is what I generally think of as the birth of linear programming as a formal discipline. Later, in 1979, Leonid Khachiyan introduced the ellipsoid method, the first algorithm proven to solve LPs in polynomial time, and in 1984 Narendra Karmarkar introduced an interior-point method that was both polynomial-time and practically fast, sparking renewed interest in LP algorithms.
Problem Statement
I define a linear program in standard form as: find a vector x that maximizes (or minimizes) a linear objective function c^T x, subject to linear inequality constraints Ax ≤ b and non-negativity constraints x ≥ 0. I want an algorithm that either finds an optimal x, reports that the problem is infeasible (no x satisfies the constraints), or reports that it is unbounded (the objective can be improved indefinitely).
Core Concepts
- Objective function: the linear function I am trying to maximize or minimize, e.g.,
c_1 x_1 + c_2 x_2 + ... + c_n x_n. - Constraints: linear inequalities or equalities that restrict feasible values of
x. - Feasible region: the set of all points satisfying every constraint, geometrically a convex polytope.
- Basic feasible solution: a vertex (corner point) of the feasible region.
- Slack variables: additional variables I introduce to convert inequality constraints into equality constraints.
- Duality: every LP (the “primal”) has an associated “dual” LP whose optimal value equals the primal’s optimal value (under mild conditions).
How It Works
I focus here on the simplex method, the classical algorithm for solving LPs:
- I convert the LP into standard form, introducing slack variables so all constraints become equalities.
- I find an initial basic feasible solution — often the origin, if that is feasible.
- I set up a “tableau” representing the current solution and the objective function in terms of the current basic variables.
- I check if any non-basic variable, if increased, would improve the objective. If none would, I have reached optimality and I stop.
- Otherwise, I choose an entering variable (one that would improve the objective) and a leaving variable (chosen via the minimum-ratio test to maintain feasibility).
- I perform a “pivot” operation, updating the tableau to reflect the new basic feasible solution.
- I repeat steps 4–6 until no improving variable exists (optimal) or I detect the problem is unbounded.
Working Principle
I think of the simplex method as walking along the edges of the feasible region’s polytope, moving from one vertex to an adjacent vertex that has an equal or better objective value, until I reach a vertex where no adjacent vertex is better. Because the feasible region is convex and the objective is linear, I know that a local optimum found this way (no improving neighbor) is also a global optimum — this is the key geometric fact that makes the simplex method correct.
Mathematical Foundation
I express the standard form of a linear program as:
$$ \text{maximize } c^T x \quad \text{subject to } Ax \le b, ; x \ge 0 $$
where $x \in \mathbb{R}^n$, $c \in \mathbb{R}^n$, $A \in \mathbb{R}^{m \times n}$, and $b \in \mathbb{R}^m$.
I convert inequality constraints to equalities by introducing a slack vector $s \ge 0$:
$$ Ax + s = b $$
The associated dual program is:
$$ \text{minimize } b^T y \quad \text{subject to } A^T y \ge c, ; y \ge 0 $$
Weak duality guarantees that for any feasible primal x and dual y:
$$ c^T x \le b^T y $$
and strong duality guarantees that at optimality:
$$ c^T x^* = b^T y^* $$
Diagrams
flowchart TD
A["Convert LP to standard form with slack variables"] --> B["Find initial basic feasible solution"]
B --> C{"Is there an improving entering variable?"}
C -- No --> D["Optimal solution found - stop"]
C -- Yes --> E["Choose entering variable"]
E --> F["Minimum ratio test to choose leaving variable"]
F --> G{"Ratio test empty?"}
G -- Yes --> H["Problem is unbounded - stop"]
G -- No --> I["Pivot: update tableau"]
I --> C
Pseudocode
SIMPLEX(A, b, c)
initialize tableau with slack variables
while there exists a non-basic variable j with reduced cost c_j > 0
choose entering variable j (e.g., most positive reduced cost)
for each row i with A[i][j] > 0
ratio[i] = b[i] / A[i][j]
if no such row exists
return "unbounded"
choose leaving variable r = row with minimum ratio
pivot(tableau, r, j) // update tableau so column j becomes a unit vector
return current basic feasible solution as optimal
Step-by-Step Example
I solve the following small LP:
$$ \text{maximize } z = 3x_1 + 5x_2 $$ $$ \text{subject to } x_1 \le 4, \quad 2x_2 \le 12, \quad 3x_1 + 2x_2 \le 18, \quad x_1, x_2 \ge 0 $$
I introduce slack variables s1, s2, s3 for the three constraints:
x1 + s1 = 42x2 + s2 = 123x1 + 2x2 + s3 = 18
Starting at the origin (x1=0, x2=0), I check reduced costs: increasing x2 improves z faster (coefficient 5 vs 3), so I choose x2 as the entering variable. The minimum-ratio test across rows gives me a leaving variable, and after pivoting I find x2 = 6. I repeat this process, and after one more pivot (bringing in x1), I reach the optimal vertex x1 = 2, x2 = 6, giving:
$$ z = 3(2) + 5(6) = 6 + 30 = 36 $$
I verify this is optimal because no further pivot improves z — every remaining reduced cost is non-positive.
Time Complexity
- Worst case: the simplex method has exponential worst-case time complexity, $O(2^n)$ in pathological constructed examples (e.g., the Klee-Minty cube), since in principle it could visit every vertex of the polytope.
- Average/practical case: despite the worst-case bound, I find the simplex method runs very efficiently in practice, often in low-polynomial time for real-world problems, which is why it remained the dominant LP algorithm for decades.
- Polynomial-time alternatives: the ellipsoid method runs in polynomial time but is slow in practice, while interior-point methods run in polynomial time, roughly $O(n^{3.5} L)$ (where
Lis the input bit-length), and are competitive with simplex on large-scale problems.
Space Complexity
I need $O(mn)$ space to store the simplex tableau, where m is the number of constraints and n is the number of variables, since the tableau essentially stores the constraint matrix A (augmented with slack variable columns) plus the objective row.
Correctness Analysis
I justify the simplex method’s correctness using convexity: the feasible region of a linear program is a convex polytope, and I can prove that a linear function attains its maximum (if bounded) at a vertex of this polytope. The simplex method only moves between adjacent vertices that do not decrease the objective, and it terminates exactly when no adjacent vertex improves the objective — at that point, because of convexity, I know I have found a vertex that is globally optimal, not just locally optimal among its neighbors. Anti-cycling rules, such as Bland’s rule, further guarantee the algorithm terminates rather than looping forever among degenerate solutions.
Advantages
- I get an algorithm that is extremely efficient in practice, even though its worst-case complexity is exponential.
- LP models are highly expressive, capable of encoding a huge range of resource-allocation, scheduling, and network problems.
- Duality theory gives me powerful tools for sensitivity analysis and for bounding solutions.
- Mature, highly optimized commercial and open-source solvers exist (e.g., CPLEX, Gurobi, GLPK), so I rarely need to implement simplex from scratch.
Disadvantages
- The simplex method has exponential worst-case time complexity, which can matter for adversarially constructed inputs.
- LP requires the objective and constraints to be strictly linear — I cannot directly model nonlinear relationships without approximation or reformulation.
- Degenerate solutions can cause cycling if I don’t apply anti-cycling rules like Bland’s rule.
- Numerical stability can be an issue with naive tableau implementations on ill-conditioned problems.
Applications
- Supply chain and logistics: minimizing transportation costs subject to supply and demand constraints.
- Manufacturing: optimizing production schedules to maximize profit given resource limits.
- Finance: portfolio optimization under risk and budget constraints.
- Telecommunications: network flow and bandwidth allocation problems.
- Agriculture: optimal crop planning subject to land, water, and labor constraints.
- Operations research more broadly: LP is often the first tool I reach for before considering more complex integer or nonlinear programming.
Implementation in C
#include <stdio.h>
#include <float.h>
#define ROWS 4 /* 3 constraints + 1 objective row */
#define COLS 7 /* x1, x2, s1, s2, s3, RHS, (placeholder) */
/* I implement a compact tabular simplex method for a small maximization LP.
Tableau layout: [x1, x2, s1, s2, s3 | RHS], last row is the objective row. */
void printTableau(double tableau[ROWS][COLS]) {
for (int i = 0; i < ROWS; i++) {
for (int j = 0; j < COLS - 1; j++) {
printf("%8.2f", tableau[i][j]);
}
printf("\n");
}
printf("\n");
}
void simplex(double tableau[ROWS][COLS]) {
int numVars = 5; /* x1, x2, s1, s2, s3 */
while (1) {
/* I look for a positive coefficient in the objective row (row 3) to improve z */
int pivotCol = -1;
double best = 1e-9;
for (int j = 0; j < numVars; j++) {
if (tableau[3][j] > best) {
best = tableau[3][j];
pivotCol = j;
}
}
if (pivotCol == -1) {
break; /* no improving column: I have reached optimality */
}
/* I run the minimum ratio test to find the pivot row */
int pivotRow = -1;
double minRatio = DBL_MAX;
for (int i = 0; i < 3; i++) {
if (tableau[i][pivotCol] > 1e-9) {
double ratio = tableau[i][5] / tableau[i][pivotCol];
if (ratio < minRatio) {
minRatio = ratio;
pivotRow = i;
}
}
}
if (pivotRow == -1) {
printf("Problem is unbounded.\n");
return;
}
/* I normalize the pivot row */
double pivotVal = tableau[pivotRow][pivotCol];
for (int j = 0; j < COLS - 1; j++) {
tableau[pivotRow][j] /= pivotVal;
}
/* I eliminate the pivot column from every other row */
for (int i = 0; i < ROWS; i++) {
if (i != pivotRow) {
double factor = tableau[i][pivotCol];
for (int j = 0; j < COLS - 1; j++) {
tableau[i][j] -= factor * tableau[pivotRow][j];
}
}
}
}
}
int main(void) {
/* I encode: maximize z = 3x1 + 5x2
subject to: x1 <= 4, 2x2 <= 12, 3x1 + 2x2 <= 18 */
double tableau[ROWS][COLS] = {
/* x1 x2 s1 s2 s3 RHS */
{ 1, 0, 1, 0, 0, 4 },
{ 0, 2, 0, 1, 0, 12 },
{ 3, 2, 0, 0, 1, 18 },
{ 3, 5, 0, 0, 0, 0 } /* objective row (I maximize, coefficients kept positive) */
};
simplex(tableau);
printf("Optimal tableau:\n");
printTableau(tableau);
printf("Optimal z = %.2f\n", -tableau[3][5] == 0 ? -tableau[3][5] : tableau[3][5]);
return 0;
}
I kept this implementation intentionally small and specific to a 2-variable, 3-constraint example so I could clearly show the tableau operations (pivot selection, ratio test, row reduction) without the complexity of a fully general simplex implementation, which would need to handle arbitrary numbers of variables and constraints, plus phase-one feasibility handling.
Sample Input and Output
Input: maximize z = 3x1 + 5x2 subject to x1 ≤ 4, 2x2 ≤ 12, 3x1 + 2x2 ≤ 18, x1, x2 ≥ 0.
Output:
Optimal tableau:
[final row-reduced tableau values]
Optimal z = 36.00which matches the hand-computed optimum I found earlier, x1 = 2, x2 = 6, z = 36.
Optimization Techniques
- Revised simplex method: I can avoid storing and updating the entire tableau by maintaining only the basis inverse, significantly reducing memory and computation for large, sparse problems.
- Bland’s rule / anti-cycling rules: I apply these to guarantee termination on degenerate problems without sacrificing much practical performance.
- Interior-point methods: for very large-scale LPs, I often switch to interior-point methods, which scale better than simplex on certain problem structures.
- Warm starting: when I re-solve a slightly modified LP (e.g., in branch-and-bound for integer programming), I reuse the previous optimal basis to speed up convergence.
- Sparse matrix representations: for large, sparse constraint matrices, I use sparse linear algebra techniques rather than dense tableau storage.
Common Mistakes
- Forgetting to convert a minimization problem correctly when using a maximization-oriented simplex implementation (I need to negate the objective, then negate the result back).
- Mishandling the ratio test — forgetting to only consider rows with a strictly positive pivot-column coefficient.
- Not detecting unboundedness, leading to infinite loops when no row qualifies for the ratio test.
- Ignoring degenerate cases that can cause cycling without an anti-cycling rule.
- Confusing the primal and dual formulations when interpreting shadow prices or sensitivity results.
Further Reading
- Cormen, Leiserson, Rivest, Stein — Introduction to Algorithms, Chapter 29 (Linear Programming): https://mitpress.mit.edu/9780262046305/introduction-to-algorithms/
- Dantzig — Linear Programming and Extensions: https://www.rand.org/pubs/reports/R366.html
- Bertsimas, Tsitsiklis — Introduction to Linear Optimization: https://www.athenasc.com/linoptbook.html
- Karmarkar — “A New Polynomial-Time Algorithm for Linear Programming,” Combinatorica, 1984: https://link.springer.com/article/10.1007/BF02579150
- Stanford CS261 — Optimization and Algorithmic Paradigms lecture notes: https://web.stanford.edu/class/cs261/
