For finitely supported measures, optimal transport is a linear program. This is not merely a numerical approximation: the finite model displays the primal polytope, dual potentials, complementary slackness, integrality, and nonuniqueness in completely explicit form.
The transport polytope
Definition. Let \[ \mu=\sum_{i=1}^m a_i\delta_{x_i}, \qquad \nu=\sum_{j=1}^n b_j\delta_{y_j}, \] where \(a_i,b_j\ge0\) and \(\sum_i a_i=\sum_j b_j=1\). A coupling is an \(m\times n\) matrix \(P=(p_{ij})\) satisfying \[ p_{ij}\ge0,\qquad \sum_jp_{ij}=a_i,\qquad \sum_ip_{ij}=b_j. \] The feasible set \[ U(a,b)=\{P\ge0:P\mathbf1_n=a,\ P^{\mathsf T}\mathbf1_m=b\} \] is the transport polytope.
Proposition. If all \(a_i,b_j>0\), the affine hull of \(U(a,b)\) has dimension \((m-1)(n-1)\). The polytope is nonempty, compact, and convex.
Proof. The product matrix \(p_{ij}=a_ib_j\) is feasible, proving nonemptiness. Closedness and boundedness in finite dimension give compactness, and the constraints are linear, giving convexity. There are \(m+n\) equality constraints, but one is redundant because total row mass equals total column mass. Their rank is \(m+n-1\), so the affine dimension is \(mn-(m+n-1)=(m-1)(n-1)\).
Primal and dual linear programs
Finite Kantorovich duality. For the cost matrix \(C=(c_{ij})\), \[ \min_{P\in U(a,b)}\sum_{i,j}c_{ij}p_{ij} = \max_{\substack{u_i+v_j\le c_{ij}}} \left\{\sum_i a_iu_i+\sum_jb_jv_j\right\}. \]
Proof. Introduce multipliers \(u_i\) and \(v_j\) for the row and column constraints. The Lagrangian is \[ \sum_{i,j}c_{ij}p_{ij} +\sum_i u_i\left(a_i-\sum_jp_{ij}\right) +\sum_j v_j\left(b_j-\sum_ip_{ij}\right) = a\cdot u+b\cdot v+\sum_{i,j}(c_{ij}-u_i-v_j)p_{ij}. \] Its infimum over \(p_{ij}\ge0\) is finite exactly when \(u_i+v_j\le c_{ij}\), in which case it equals \(a\cdot u+b\cdot v\). Finite-dimensional linear-programming duality gives equality of the two optimal values because the primal polytope is nonempty and bounded.
Complementary slackness. A feasible \(P\) and feasible \((u,v)\) are optimal if and only if \[ p_{ij}>0\quad\Longrightarrow\quad u_i+v_j=c_{ij}. \] Thus positive transport occurs only along tight dual edges.
Assignment and the Birkhoff polytope
Birkhoff--von Neumann theorem. Every \(n\times n\) doubly stochastic matrix is a convex combination of permutation matrices. Consequently, the extreme points of \[ \mathcal B_n=\{P\ge0:P\mathbf1=\mathbf1,\ P^{\mathsf T}\mathbf1=\mathbf1\} \] are exactly the permutation matrices.
Proof. Build a bipartite graph with edge \((i,j)\) whenever \(p_{ij}>0\). For any set \(S\) of row vertices, the total mass leaving \(S\) equals \(|S|\), while each neighboring column receives at most one; hence Hall's condition \(|N(S)|\ge|S|\) holds. There is a perfect matching \(\sigma\) supported on positive entries. Let \(\varepsilon=\min_i p_{i,\sigma(i)}>0\). Then \[ P=\varepsilon P_\sigma+(1-\varepsilon)Q, \] where \(Q=(P-\varepsilon P_\sigma)/(1-\varepsilon)\) is again doubly stochastic. At least one additional entry becomes zero. Iteration terminates after finitely many steps and gives a convex combination of permutation matrices. A permutation matrix cannot be written as a nontrivial convex combination of other doubly stochastic matrices because its zero pattern forces both matrices in the combination to coincide with it.
Assignment corollary. For equal-weight measures with \(m=n\), a linear transport cost has an optimal permutation plan. The Kantorovich relaxation therefore has the same value as the discrete assignment problem, though the relaxed optimizer may be nonunique.
Network structure and integrality
The constraint matrix is the incidence matrix of a bipartite graph and is totally unimodular. Hence if the supplies and demands are integers, every vertex of the transport polytope is integral. This is why a min-cost flow algorithm returns an integral shipment without adding integer constraints. In the equal-unit case, integrality reduces the vertices to assignments.
Sparse vertices. Every nondegenerate vertex of \(U(a,b)\) has exactly \(m+n-1\) positive entries; every vertex has at most \(m+n-1\) positive entries.
Proof. Associate a bipartite support graph to \(P\). If its positive edges contain a cycle, add and subtract a small \(\varepsilon\) alternately around the cycle. Both resulting matrices remain feasible and distinct, with \(P\) their midpoint, so \(P\) is not extreme. Thus a vertex has acyclic support. A forest on \(m+n\) vertices has at most \(m+n-1\) edges. Nondegeneracy means the support is a spanning tree, giving equality.
A complete two-by-three example
Example. Take \[ a=(1/2,1/2),\qquad b=(1/4,1/2,1/4),\qquad C=\begin{pmatrix}0&1&3\\2&1&0\end{pmatrix}. \] Consider \[ P^\star= \begin{pmatrix} 1/4&1/4&0\\ 0&1/4&1/4 \end{pmatrix}. \] Its row and column sums are \(a,b\), and its cost is \(1/2\). The dual variables \[ u=(0,0),\qquad v=(0,1,0) \] satisfy \(u_i+v_j\le c_{ij}\) and have dual value \(b\cdot v=1/2\). Primal and dual values agree, proving optimality. Every positive entry of \(P^\star\) lies on a tight edge, which displays complementary slackness directly.
From discrete to continuous transport
Approximating \(\mu,\nu\) by empirical or quantized measures turns continuous transport into a sequence of finite linear programs. The stability theorem from page 3 explains when their optimal values and plans converge. Entropic regularization will later replace the sparse polytope optimizer by a smooth positive matrix that can be found through diagonal scaling.