Min-Cost Flow via Cost Scaling, in depth: epsilon-optimality, the refine phase, a tested implementation and the heuristics that make it fast

Where cost scaling fits

Successive shortest paths may need as many augmentations as the total supply. Network simplex is often fastest on sparse, moderate instances but has no clean polynomial bound. Cost scaling's running time depends on the logarithm of the largest cost C, not on the flow value. The intuition is that of any scaling algorithm: a coarse tolerance is cheap to reach, and tightening it by a constant factor needs only local repairs because the previous solution is nearly right. After log(nC) phases integer costs leave no room for error, and you stop with an exact optimum.

Prices, reduced costs and epsilon-optimality

Give every node v a price p(v). The reduced cost of a residual arc (u, v) with cost c is c_p(u,v) = c(u,v) + p(u) - p(v). Around any cycle the prices cancel, so a cycle's reduced cost equals its true cost. A flow is optimal exactly when no residual cycle has negative cost, which is equivalent to the existence of prices with every residual reduced cost at least zero (the complementary slackness conditions of the LP).

Cost scaling relaxes this. A flow is epsilon-optimal for prices p if every residual arc has c_p >= -eps. Any flow is epsilon-optimal for a large enough epsilon: with zero prices, eps = C (the largest absolute cost) works. The key lemma tells you when to stop. A simple residual cycle has at most n arcs, so its cost is at least -n * eps. If costs are integers and eps < 1/n, the cycle cost is greater than -1 and therefore at least 0: no negative cycle exists, and the flow is optimal.

Fractional epsilon is awkward, so implementations multiply every cost by alpha = n + 1 at the start. On the scaled costs, eps = 1 means eps < 1/n on the original costs, and every quantity stays an integer. The outer loop is then three lines: set eps to the largest scaled cost, repeatedly divide eps by a factor (2 in the code below, 8 to 16 in tuned libraries) and call refine, and stop after the phase with eps = 1.

Refine: repairing a pseudoflow with push and relabel

Refine receives a flow that is 2eps-optimal (or alpha-eps-optimal) and must return one that is eps-optimal. It does this in two moves. First, it saturates every residual arc whose reduced cost is negative. After that, every residual arc has c_p at least 0, so the flow is trivially 0-optimal, but it is no longer a flow: conservation is broken and some nodes have excess (more in than out) while others have deficit. It is a pseudoflow.

Second, it repairs conservation with push and relabel, exactly as push-relabel max flow does, but with reduced costs playing the role of heights. An arc is admissible if it has residual capacity and c_p < 0. A node with positive excess pushes along admissible arcs. When it has excess and no admissible arc, it is relabelled: its price drops to max over residual (u,w) of p(w) - c(u,w) - eps, which makes the best outgoing arc's reduced cost exactly -eps. Both operations preserve eps-optimality: a push creates the reverse arc with reduced cost greater than 0, and the relabel is chosen so no arc out of u falls below -eps (arcs into u only get more expensive in reduced terms, which is safe).

Why does this terminate? Pushes only travel downhill in reduced cost, so they cannot cycle without a relabel, and Goldberg and Tarjan showed that in one refine each price can fall by at most O(n eps) in total, so each node is relabelled O(n) times. That gives O(n^2) relabels and, for the generic method, O(n^2 m) pushes per phase. The bound depends on a feasible instance: if excess cannot reach any deficit, refine never ends, which is why the wrapper below adds an expensive bypass arc.

A tested implementation

The implementation stores each arc next to its reverse (arc e and e ^ 1), keeps a FIFO queue of active nodes and a current-arc pointer per node, and counts operations so you can see the work each phase does. The wrapper routes a demand from s to t; it adds a bypass arc s to t whose cost exceeds the cost of any simple path, so the instance is always feasible and the bypass carries flow only if the real network cannot.

from collections import deque

class CostScalingMCF:
    def __init__(self, n):
        self.n = n
        self.to, self.cap, self.cost, self.adj = [], [], [], [[] for _ in range(n)]

    def add_edge(self, u, v, cap, cost):          # arc e and its reverse e ^ 1
        for a, b, c, w in ((u, v, cap, cost), (v, u, 0, -cost)):
            self.adj[a].append(len(self.to))
            self.to.append(b); self.cap.append(c); self.cost.append(w)
        return len(self.to) - 2

    def solve(self, supply):
        alpha = self.n + 1                         # eps < 1 on scaled costs => optimal
        cost = [w * alpha for w in self.cost]
        excess, price = list(supply), [0] * self.n
        eps = max([abs(w) for w in cost] + [1])
        stats = {"phases": 0, "pushes": 0, "relabels": 0}
        while True:
            eps = max(1, eps // 2)
            stats["phases"] += 1
            self._refine(eps, cost, excess, price, stats)
            if eps == 1:
                break
        total = sum(self.cost[e] * self.cap[e ^ 1] for e in range(0, len(self.to), 2))
        return total, price, stats

    def _refine(self, eps, cost, excess, price, stats):
        to, cap, adj = self.to, self.cap, self.adj
        for u in range(self.n):                    # saturate negative reduced costs
            for e in adj[u]:
                v = to[e]
                if cap[e] > 0 and cost[e] + price[u] - price[v] < 0:
                    f = cap[e]
                    cap[e] -= f; cap[e ^ 1] += f
                    excess[u] -= f; excess[v] += f
        active = deque(u for u in range(self.n) if excess[u] > 0)
        current = [0] * self.n
        while active:
            u = active.popleft()
            while excess[u] > 0:
                if current[u] == len(adj[u]):      # no admissible arc left: relabel
                    price[u] = max(price[to[e]] - cost[e] for e in adj[u] if cap[e] > 0) - eps
                    current[u] = 0
                    stats["relabels"] += 1
                    continue
                e = adj[u][current[u]]
                v = to[e]
                if cap[e] > 0 and cost[e] + price[u] - price[v] < 0:   # admissible: push
                    f = min(excess[u], cap[e])
                    cap[e] -= f; cap[e ^ 1] += f
                    excess[u] -= f; excess[v] += f
                    stats["pushes"] += 1
                    if 0 < excess[v] <= f:         # v just became active
                        active.append(v)
                else:
                    current[u] += 1

def min_cost_flow(n, edges, s, t, demand):
    big = n * max([abs(w) for *_, w in edges] + [1]) + 1
    g = CostScalingMCF(n)
    for u, v, cap, cost in edges:
        g.add_edge(u, v, cap, cost)
    bypass = g.add_edge(s, t, demand, big)
    supply = [0] * n; supply[s] = demand; supply[t] = -demand
    total, _, stats = g.solve(supply)
    if g.cap[bypass ^ 1] > 0:
        return None, stats                         # infeasible: the bypass carried flow
    return total, stats

The relabel takes the maximum only over arcs with residual capacity; a node with excess always has one, because the excess arrived along an arc whose reverse is now residual.

Worked example: four nodes, three demands

Take four nodes, s = 0 and t = 3, with arcs written as (from, to, capacity, cost): (0,1,3,2), (0,2,2,5), (1,2,2,1), (1,3,2,6), (2,3,3,2). Ask for 4 units.

Check the answer by hand first. The cheapest path is 0-1-2-3 at cost 5, limited to 2 units by arc (1,2). Arc (2,3) then has 1 unit left, so 0-2-3 at cost 7 carries 1 unit. The last unit goes 0-1-3 at cost 8. Total 2 * 5 + 7 + 8 = 25. The code returns (25, {'phases': 6, 'pushes': 30, 'relabels': 19}): with n + 1 = 5 the largest scaled cost is the bypass (4 * 6 + 1 = 25, times 5 = 125), so eps runs 62, 31, 15, 7, 3, 1, six phases.

Ask for 5 units and it returns 35: the fifth unit takes 0-2, the reverse of (1,2) at cost -1, then 1-3, for 10 more, a reroute that appears here as a push on a reverse arc. Ask for 6 and it returns None, because only 5 units can leave s; the bypass absorbed the sixth, which is how the wrapper detects infeasibility without a separate max-flow run.

Against a Bellman-Ford SSP on 3,000 random graphs (2 to 8 nodes, up to 16 arcs), the code gave zero mismatches, infeasible cases included. On a random 400-node, 6,000-arc instance with demand 200 the Python versions took 0.74 s and 0.80 s, a tie; the advantage is asymptotic, in large supplies, not in small-graph timings.

Complexity and the heuristics that matter

With the generic selection rule, one refine costs O(n^2 m) and there are O(log(nC)) phases, so the whole algorithm is O(n^2 m log(nC)). Processing active nodes in FIFO or wave order brings a phase to O(n^3), and Goldberg and Tarjan's dynamic-tree version reaches O(nm log(n^2/m) log(nC)). None of these depend on the total supply, which is the point.

Real implementations are fast because of heuristics rather than the bounds:

• A larger scaling factor. Dividing eps by 8 to 16 instead of 2 means fewer phases at the price of more work per phase; LEMON's CostScaling defaults to 16.
• Price updates. The analogue of the global relabel in max flow: periodically run a Dijkstra-like pass from the deficit nodes to raise prices in one sweep instead of many small relabels.
• Arc fixing. Arcs whose reduced cost is far above eps can be dropped from scans.
• Push look-ahead and partial augmentation. Push only if the target can pass the flow on, and follow short admissible paths instead of single arcs (LEMON's PARTIAL_AUGMENT).

Operational guidance

• Use a library. LEMON (C++) ships CostScaling alongside NetworkSimplex and CapacityScaling; OR-Tools' min_cost_flow.h is a cost-scaling push-relabel solver. Write your own only to learn.
• Watch integer width. Costs are multiplied by about n, and prices can drift to around n times the largest scaled cost. With 64-bit integers and n near a million, costs above roughly 10^6 can overflow; check the product before you solve.
• Integerise costs deliberately. The termination lemma needs integer costs. Scale money to cents and round consistently; never feed floats.
• Keep the prices. Final prices are a dual certificate you can verify, and a warm start when only a few arcs change.

Failure modes

• Infinite loop on infeasible input. Refine assumes excess can reach a deficit. Without a bypass arc or a prior max-flow check, an infeasible instance spins forever.
• Wrong stopping rule. Stopping at eps = 1 on unscaled costs leaves the flow 1-optimal, not optimal; the multiplier n + 1 is what makes eps = 1 sufficient.
• Relabelling with saturated arcs. Taking the max over all arcs, including those with zero residual capacity, gives a price that can make a non-residual arc look admissible and corrupts the invariant.
• Missing re-activation. Forgetting to enqueue a node whose excess turns positive leaves excess stranded and the result is not a flow; assert conservation at the end.

Trade-offs

What to do next

1. Run the code on the worked example, print prices after each phase, and confirm every residual reduced cost is at least -eps at the end of each refine.
2. Add an assertion that excess is zero everywhere after the last phase, and a randomized comparison against an SSP reference such as the one in the min-cost flow article.
3. Change the scaling factor from 2 to 8 and 16 and record phases, pushes and relabels on a few thousand-node instances.
4. Model one real assignment problem as min-cost flow and solve it both ways; compare against the Hungarian algorithm on dense cases.
5. Revisit Bellman-Ford and its potentials: the prices here are the same dual objects, maintained incrementally.
6. Before production use, try LEMON's CostScaling or OR-Tools on your data and verify the dual certificate rather than trusting the solver's status flag alone.