Multi-Commodity Flow, in depth: LP formulations, integrality, the flow-cut gap and Garg-Konemann

The problem and its three objectives

Take a directed graph G = (V, E) with capacities c(e) ≥ 0, and k commodities. Commodity i has a source s_i, a sink t_i and, for some objectives, a demand d_i. A multi-commodity flow assigns each commodity its own flow f_i(e) ≥ 0 such that:

• conservation per commodity: at every vertex other than s_i and t_i, commodity i's inflow equals its outflow;
• joint capacity: on every edge, the sum over commodities, Σ_i f_i(e), is at most c(e).

Three objectives come up repeatedly. Maximum total flow maximises the sum of all commodities' values and does not care who gets what. Maximum concurrent flow maximises the largest λ such that every commodity ships at least λ·d_i; it is the fairness objective, and λ < 1 means the network cannot meet all demands. Minimum-cost multi-commodity flow ships every demand at minimum total cost, generalising single-commodity min-cost flow.

The crucial point: you cannot solve this by running k independent max-flows. Each would happily use the full capacity of a shared edge, and their sum would overflow it. Nor can you add a super-source and super-sink as in multi-terminal single-commodity flow, because that lets flow from s₁ end at t₂.

A worked example by hand

Work the example by hand before reaching for a solver. Commodity 1 can use the private path s1-c-t1 (capacity 1) and the shared path s1-a-b-t1. Commodity 2 can use only s2-a-b-t2. Edge a-b has capacity 2.

• Each commodity alone: commodity 1 could ship 1 + 2 = 3 and commodity 2 could ship 2. The naive sum of 5 is infeasible.
• Maximum total: the private path carries 1, and a-b carries 2 of anything, so the total is 3. Many splits achieve it, including one that starves commodity 2 completely.
• Maximum concurrent with demands 2 and 2: let x be commodity 1's share of a-b and y commodity 2's share, with x + y ≤ 2. We need 1 + x ≥ 2λ and y ≥ 2λ. Balancing gives x = 0.5, y = 1.5, λ = 0.75: both commodities get 1.5 units, 75% of their demand.

Note that the fair answer gives commodity 1 a fractional half-unit on a-b. That is not an accident of this example, and it matters a great deal below.

The LP formulations, with code

The edge formulation has one variable per commodity per edge, k·|E| in total, one conservation equality per commodity per internal vertex, and one joint capacity inequality per edge. Any LP solver handles it directly. The code below builds it for SciPy's linprog with the HiGHS backend and supports both the total and the concurrent objective.

import numpy as np
from scipy.optimize import linprog

def mcf_lp(nodes, edges, commodities, concurrent=False):
    # Edge formulation. edges: list of (u, v, cap); commodities: list of (s, t, demand).
    m, k = len(edges), len(commodities)
    nv = m * k + (1 if concurrent else 0)          # f[i,e] laid out commodity-major, then lambda
    idx = lambda i, e: i * m + e
    c = np.zeros(nv)
    A_ub, b_ub, A_eq, b_eq = [], [], [], []
    for e, (_, _, cap) in enumerate(edges):        # joint capacity: sum_i f[i,e] <= cap
        row = np.zeros(nv)
        for i in range(k):
            row[idx(i, e)] = 1
        A_ub.append(row); b_ub.append(cap)
    for i, (s, t, d) in enumerate(commodities):
        for v in nodes:                            # conservation everywhere except s and t
            if v in (s, t):
                continue
            row = np.zeros(nv)
            for e, (a, b, _) in enumerate(edges):
                if a == v: row[idx(i, e)] += 1
                if b == v: row[idx(i, e)] -= 1
            A_eq.append(row); b_eq.append(0)
        out = np.zeros(nv)                         # net outflow of s = value of commodity i
        for e, (a, b, _) in enumerate(edges):
            if a == s: out[idx(i, e)] += 1
            if b == s: out[idx(i, e)] -= 1
        if concurrent:                             # value_i >= lambda * d_i
            row = -out; row[-1] = d
            A_ub.append(row); b_ub.append(0)
        else:
            c -= out                               # maximise total value
    if concurrent:
        c[-1] = -1
    res = linprog(c, A_ub=np.array(A_ub), b_ub=b_ub,
                  A_eq=np.array(A_eq) if A_eq else None, b_eq=b_eq or None,
                  bounds=[(0, None)] * nv, method="highs")
    return -res.fun, res.x

nodes = ["s1", "s2", "a", "b", "c", "t1", "t2"]
edges = [("s1", "a", 2), ("s2", "a", 2), ("a", "b", 2), ("b", "t1", 2),
         ("b", "t2", 2), ("s1", "c", 1), ("c", "t1", 1)]
comms = [("s1", "t1", 2), ("s2", "t2", 2)]
print(mcf_lp(nodes, edges, comms)[0])                    # 3.0
print(mcf_lp(nodes, edges, comms, concurrent=True)[0])   # 0.75

Run against the example, it prints a maximum total of 3.0 and λ = 0.75, and the concurrent solution routes 0.5 of commodity 1 and 1.5 of commodity 2 over a-b, exactly the hand calculation. Dense rows are fine for teaching; for real graphs build scipy.sparse matrices, because the constraint matrix has about k·(|V| + |E|) rows.

The path formulation instead has one variable per commodity per s-t path. It has exponentially many columns, but optimal solutions use few of them, so it is solved by column generation: solve a restricted LP, read the dual price y(e) of each capacity constraint, and for each commodity run a shortest-path search under edge lengths y(e). A path shorter than the commodity's dual threshold improves the LP; add it and repeat. The pricing step is just Dijkstra, which is why this is the formulation production solvers use. The duality behind it is covered in LP duality.

What breaks: integrality and the flow-cut gap

Single-commodity flow has two beautiful properties: integer capacities give an integer optimal flow, and the maximum flow equals the minimum cut. Both fail here.

• Integrality. Finding an integral multi-commodity flow, for example routing each demand as whole units, is NP-hard even with two commodities (Even, Itai and Shamir, 1976). The fractional LP is polynomial, but its optimum may be fractional, as the 0.5 above shows. If you need integers, you need rounding, branch-and-bound or a heuristic, and you should expect a gap.
• Max-flow min-cut. The natural analogue of a cut for concurrent flow is the sparsest cut: the capacity crossing a cut divided by the demand it separates. Concurrent flow is always at most the sparsest-cut ratio, but they can differ. Leighton and Rao showed the gap is O(log n) for uniform demands, and later work (Linial, London and Rabinovich; Aumann and Rabani) extended this to O(log k) for k arbitrary commodities. The bound is tight: on constant-degree expanders the gap really is logarithmic.
• One positive exception. For two commodities in an undirected graph with integer capacities, Hu's theorem (1963) restores max-flow min-cut, with a half-integral optimal flow.

The practical consequence is that a multi-commodity solver is an LP solver, and the cut-based reasoning you use for single-commodity bottleneck analysis only gives bounds. That flow-cut gap is also why multi-commodity flow is the engine behind approximation algorithms for graph partitioning.

Garg-Konemann: approximation with shortest paths

When the graph has millions of edges and thousands of commodities, even a sparse LP becomes heavy. Garg and Könemann's algorithm, a multiplicative-weights method, computes a (1 − ε)-approximate maximum total flow using only shortest-path computations. The idea: give every edge a length that grows exponentially with its load, repeatedly push flow along the globally shortest source-sink path, and stop when even that path has become "long". Congested edges get expensive, so later flow routes around them.

import heapq, math

def garg_konemann(edges, commodities, eps=0.1):
    # Approximate maximum total multicommodity flow with multiplicative weights.
    m = len(edges)
    delta = (1 + eps) * ((1 + eps) * m) ** (-1 / eps)
    length = [delta / cap for (_, _, cap) in edges]
    adj = {}
    for e, (u, v, _) in enumerate(edges):
        adj.setdefault(u, []).append((e, v))

    def shortest(s, t):
        dist, prev, pq = {s: 0.0}, {}, [(0.0, s)]
        while pq:
            d, u = heapq.heappop(pq)
            if u == t: break
            if d > dist.get(u, math.inf): continue
            for e, v in adj.get(u, []):
                if d + length[e] < dist.get(v, math.inf):
                    dist[v], prev[v] = d + length[e], e
                    heapq.heappush(pq, (dist[v], v))
        if t not in dist: return math.inf, []
        path, v = [], t
        while v != s:
            path.append(prev[v]); v = edges[prev[v]][0]
        return dist[t], path

    flow, value = [0.0] * m, 0.0
    while True:
        d, path = min((shortest(s, t) for s, t, _ in commodities), key=lambda r: r[0])
        if d >= 1: break                                   # every path is now "long"
        u = min(edges[e][2] for e in path)                 # bottleneck capacity
        for e in path:
            flow[e] += u
            length[e] *= 1 + eps * u / edges[e][2]
        value += u
    scale = max(f / cap for f, (_, _, cap) in zip(flow, edges))   # make it feasible
    return value / scale, [f / scale for f in flow]

The raw flow overshoots capacities by a logarithmic factor. Dividing by the observed maximum congestion is always feasible; the published analysis fixes its own scaling constant, so treat this code as a practical variant and check it against the LP, as here. On the worked example it returned 2.8333 at ε = 0.3, 2.9321 at ε = 0.1 and 2.9663 at ε = 0.05, converging on the LP optimum of 3 from below, with every edge at or under capacity. The number of iterations grows roughly as (m/ε²)·log m, so tightening ε is the main cost knob. Variants of the same scheme handle concurrent flow and min-cost flow.

Using it in production

Multi-commodity flow sits at the core of wide-area traffic engineering: systems such as Microsoft's SWAN and Google's B4 allocate bandwidth between site pairs with multi-commodity formulations, usually with fairness layers on top and approximate solvers for speed. It also appears in freight and airline network design, VLSI global routing and the analysis of interconnects in large training clusters, where commodities are collective-communication flows between accelerators. A few habits carry across all of them:

• Aggregate commodities by source. Commodities that share a source can be merged into one flow with several sinks when you only need edge loads, cutting variables from k·|E| to |V|·|E| or less.
• Restrict paths. Production traffic engineering usually precomputes a handful of candidate tunnels per pair and solves the path formulation over them. You lose a little optimality and gain predictable solve times and deployable routes.
• Pick the objective deliberately. Maximum total flow starves commodities; maximum concurrent flow is fair but rigid; max-min fairness, solved as a sequence of concurrent LPs, is what operators typically want.
• Warm-start. Demands change gradually. Reusing the previous basis or path set turns minutes into seconds.

Failure modes

• Treating fractional flow as routable. 0.5 units on a path may be fine for bandwidth split across ECMP hashing, and meaningless for a truck. Round, then re-check capacities.
• Summing single-commodity results. The classic bug: independent max flows overbook shared links.
• Using a super-source. It lets commodities swap destinations, silently giving a too-optimistic answer.
• Reading a min cut as the bottleneck. With several commodities the cut bound can be loose by a logarithmic factor; the LP dual prices are the real congestion signal.
• Numerical tolerance. Solvers return 1e-9 flows on unused edges. Threshold before extracting paths, or your path decomposition gains phantom routes.

Trade-offs

What to do next

1. Run mcf_lp on the example and change the capacity of a-b to 3; predict the new λ by hand before checking (it should rise to 1).
2. Swap the dense rows for scipy.sparse.lil_matrix and time it on a random graph with 1,000 nodes and 20 commodities.
3. Implement path extraction: decompose each commodity's edge flow into paths and check that the path flows sum to the commodity value.
4. Compare Garg-Könemann against the LP at several ε values on the same graph and plot error against run time.
5. Read about the dual: Lagrangian relaxation applied to the joint capacity constraints decomposes the problem into k independent shortest-path or min-cost-flow problems.