Maximum Flow and Minimum Cut

All the functions described in this section apply to general graphs whose edges are given a capacity, that is, networks. If the graph under consideration is not capacitated, then all its edges are assumed to have capacity one. To assign capacities to the edges of a graph, see Subsection Edge Decorations. To create a Magma network object (with type GrphNet) see Section Construction of Networks.

The fundamental network flow problem is the minimum cost flow problem, that is, determining a maximum flow at minimum cost from a specified source to a specified sink. Specializations of this problem are the shortest path problem (where there is no capacity constraint) which is covered in Section Distances, Shortest Paths and Minimum Weight Trees, and the maximum flow problem (where there is no cost constraint). Some of the related problems are the minimum spanning tree problem, the matching problem (a pairing of the edges of the graph according to some criteria), and the multicommodity flow problem (where arcs may carry several flows of different nature). For a comprehensive monograph on network problems, their implementation and applications, see [AMO93].

Let G be a general graph or multigraph whose edges are capacitated. If G is undirected then consider any edge {u, v} with capacity c as being equivalent to two directed edges [u, v] and [v, u], both with capacity c. We may thus assume without loss of generality that G is directed, and from now on G is called a network. Let V and E be G's vertex-set and edge-set respectively, and for every edge [u, v] in G, denote its capacity by c(u, v). If there is no edge [u, v] in E we assume that c(u, v) = 0. By convention the capacity of an edge is a non-negative integer.

Distinguish two vertices of G, the source s and the sink t. A flow in G is an integer-valued function f : V x V -> Z that satisfies the following properties:

(i)
Capacity constraint: forall u, v ∈V, f(u, v) ≤c(u, v).

(ii)
Skew symmetry: forall u, v ∈V, f(u, v) = - f(v, u).

(iii)
Flow conservation: forall u ∈V - {s, t}, ∑v∈V f(u, v) = 0.

Note that the flow could have been defined as a real-valued function if the edge capacity had been defined as real-valued.

The quantity f(u, v), which can be positive or negative, is called the net flow from vertex u to vertex v. The value of a flow f is defined as F = ∑v∈V f(s, v), that is, the total net flow out of the source.

In the maximum flow problem we wish to find a flow of maximum value from s to t. A cut of the network G with source s and sink t is a partition of V into S and T such that s ∈S and t∈T. The capacity of the cut c(S) determined by S is the sum of the capacity c(u, v) of the edges [u, v] such that u ∈S and v ∈T. It is easy to see that F ≤c(S), that is, the value of the flow from s into t cannot be larger than the capacity of any cut defined by s and t. In fact, F = c(S) if and only if F is a maximum flow and S is a cut of minimum capacity.

We have implemented two algorithms for finding the maximum flow in a network G from a source s to a sink t. One is the Dinic algorithm, the other is an instance of the generic push-relabel algorithm. We have already encountered both algorithms in Subsections Maximum Matching in Bipartite Graphs and General Vertex and Edge Connectivity in Graphs and Digraphs.

The Dinic algorithm

The Dinic algorithm consists of two phases. In the first phase, one constructs a layered network which consists of all the "useful" edges of G: A useful edge [u, v] has the property that f(u, v) < c(u, v). The first and the last layer of this layered network always contain s and t respectively. If such a layered network cannot be constructed then the flow in G is maximum and the algorithm ends.

In the second phase of the Dinic algorithm, one finds a maximal flow by constructing paths from s to t (using a depth-first search) in the layered network. Note that a maximal flow may not be maximum: A maximal flow satisfies the condition that for every path P from s to t in the layered network, there is at least a saturated edge in P. A saturated edge is one for which its flow equals its capacity. For more details on the Dinic algorithm, see [Eve79].

The Dinic algorithm has a theoretical complexity of O(|V|2 |E|), which can be improved to O(|E|3/2) if G is a zero-one network (the capacities of the edges are either 0 or 1), or O(|V|2/3 |E|) if G is a zero-one network with no parallel edges. Prior to the advent of the push-relabel method (which we describe below), the Dinic algorithm was shown to be superior in practice to other methods. The magma implementation of the Dinic algorithm only outperforms the push-relabel method for very sparse networks (a small number of edges relative to the number of vertices) whose edges have a small capacity and where the maximum flow is small. In other words, the Dinic method performs best on zero-one and very sparse networks.

The push-relabel method

The generic push-relabel algorithm does not construct a flow by constructing paths from s to v. Rather, it starts by pushing the maximum possible flow out from the source s into the neighbours of s, then pushing the excess flow at those vertices into their own neighbours. This is repeated until all vertices of G except s and t have an excess flow of zero (that is, the flow conservation property is satisfied). Of course this might mean that some flow is pushed back into s.

At initialisation all vertices are given a height of 0 except s which will have height |V|, and flow is pushed out of s. The maximum flow that can be pushed out of s is ∑v ∈V c(s, v). The idea is to push flow into vertices only in a "downward" manner, that is, from a vertex with height h + 1 into a vertex with height h. If a vertex in V - {s, t} has height h, it will be relabelled if it has an excess flow and none of its neighbours to which some flow could be pushed has height h - 1. Its new label (height) will be l + 1 where l is the height of the lowest labelled neighbour towards which flow can be pushed. Flow may only be pushed along an edge [u, v] while the capacity constraint f(u, v) ≤c(u, v) is satisfied.

So the general notion is of pushing flow "downward" and "discharging" the vertices in V - {s, t} from all non-zero excess flow they may have accumulated, a process which may require "lifting" some vertices so that discharging can proceed. It is easy to show that the height of any vertex in V - {s, t} is bounded by 2|V| - 1, so the algorithm must terminate. When this happens, the flow is at a maximum.

The theoretical complexity of the generic push-relabel algorithm is O(|V|2 |E|) which can be improved to at least O(|V||E| log(|V|2/|E|) thanks to some heuristics.

One commonly used heuristic is the global relabelling procedure whereby the vertices in V - {s, t} are pre-labelled according to their distance from the sink. Global relabelling may occur at specified intervals during execution. Another very useful heuristic is gap relabelling which involves recognising those vertices that have become unreachable from the sink and labelling them accordingly so that they won't be chosen during later execution. Finally, another heuristic involves the order in which the vertices to be discharged are processed.

When G is a one-zero network, choosing the vertex with smallest height works best, while choosing the vertex with largest height is recommended for general networks. These heuristics are fully described and analysed by B.V. Cherkassky and A.V. Goldberg et al. in [CGM+98] and [CG97]. The magma implementation incorporates all the above heuristics together with some new ones.

As a rule, whenever the push-relabel algorithm is applied to a zero-one network (such as occurs when computing the connectivity of a graph in General Vertex and Edge Connectivity in Graphs and Digraphs or finding a maximum matching in a bipartite network in Maximum Matching in Bipartite Graphs) the vertex with smallest height is chosen as the next vertex to be processed. When the algorithm is used to compute a maximum flow in a network (as below) the vertex with largest height is chosen.

For the two functions MinimumCut and MaximumFlow described below, the PushRelabel algorithm is used. It is only in the case of a very sparse zero-one network that Dinic may outperform PushRelabel.

MinimumCut(s, t : parameters) : GrphVert, GrphVert -> SeqEnum, RngIntElt
    Al: MonStgElt                       Default: "PushRelabel"
Given vertices s and t of a network G, this function returns the subset S defining a minimum cut { S, T} of V, where s ∈S, t ∈T, corresponding to the maximum flow F from s to t. The subset S is returned as a sequence of vertices of G. The maximum flow F is returned as a second value. The parameter Al enables the user to select the algorithm to be used: Al := "PushRelabel" or Al := "Dinic".
MinimumCut(Ss, Ts : parameters) : [ GrphVert ], [ GrphVert ] -> SeqEnum, RngIntElt
    Al: MonStgElt                       Default: "PushRelabel"
Given two sequences Ss and Ts of vertices of a network G, this function returns the subset S defining a minimum cut { S, T} of V, such that Ss ⊆S and Ts ⊆T, and which corresponds to the maximum flow F from the vertices in Ss to the vertices in Ts. The subset S is returned as a sequence of vertices of G. The maximum flow F is returned as a second value. The parameter Al enables the user to select the algorithm to be used: Al := "PushRelabel" or Al := "Dinic".
MaximumFlow(s, t : parameters) : GrphVert, GrphVert -> RngIntElt, SeqEnum
    Al: MonStgElt                       Default: "PushRelabel"
Given vertices s and t of a network G, this function returns the maximum flow F from s to t. The subset S defining a minimum cut { S, T} of VertexSet(N), s ∈S, t ∈T, corresponding to F is returned as the second value. The subset S is returned as a sequence of vertices of G. The parameter Al enables the user to select the algorithm to be used: Al := "PushRelabel" or Al := "Dinic".
MaximumFlow(Ss, Ts : parameters) : [ GrphVert ], [ GrphVert ] -> RngIntElt, SeqEnum
    Al: MonStgElt                       Default: "PushRelabel"
Given two sequences Ss and Ts of vertices of a network G, this function returns the maximum flow F from the vertices in Ss to the vertices in Ts. The subset S defining a minimum cut { S, T} of VertexSet(N), Ss ⊆S and Ts ⊆T, corresponding to F is returned as the second value. The subset S is returned as a sequence of vertices of G. The parameter Al enables the user to select the algorithm to be used: Al := "PushRelabel" or Al := "Dinic".
Flow(e) : GrphEdge -> RngIntElt
Given an edge e in a network G, this function returns the flow of the edge e as an integer. In this instance we require that the edges of G be explicitly assigned capacities (see Subsection Edge Decorations or Section Construction of Networks). Edges of G will carry a flow only if a flow has been constructed from a source to a sink in G. If no such flow has yet been constructed, all edges will have zero flow.
Flow(u, v) : GrphVert, GrphVert -> RngIntElt
Let u and v be adjacent vertices of a network G whose edges have been explicitly assigned capacities (see Subsection Edge Decorations or Section Construction of Networks). Flow(u, v) returns the total net flow from u to v as an integer. The total net flow from u to v is defined as the total outgoing flow from u into v minus the total ingoing flow into u from v. Thus the flow satisfies the skew symmetry Flow(u, v) = - Flow(v, u).

Edges of G will carry a flow only if a flow has been constructed from a source to a sink in G. If no such flow has yet been constructed, all edges will have zero flow.

Example Network_Flow (H160E4)

We illustrate the maximum flow algorithm by applying it to

find a maximum matching in a bipartite graph. This example exactly replicates the implementation of the maximum matching algorithm MaximumMatching).

We construct a bipartite graph.

> G := Graph< 7 | [ {5, 7}, {6, 7}, {4, 5, 6},
> {3}, {1, 3}, {2, 3}, {1, 2} ] : SparseRep := true >;
> assert IsBipartite(G);
> P := Bipartition(G);
> P;
[
  { 1, 2, 3 },
  { 4, 5, 6, 7 }
]
We add two extra vertices to G, the source s and the sink t and we construct three sets of pairs of vertices. The first contains all the pairs [s, u] where u∈P[1], the second contains all the pairs [u, t] where u∈P[2], and the third contains all the pairs [u, v] where u∈P[1] and v∈P[2]. From these sets we construct a multidigraph with capacitated edges.
> G +:= 2;
> G;
Graph
Vertex  Neighbours
1       7 5 ;
2       7 6 ;
3       6 5 4 ;
4       3 ;
5       3 1 ;
6       3 2 ;
7       2 1 ;
8       ;
9       ;
> s := Order(G)-1;
> t := Order(G);
>
> E1 := { [s, Index(u)] : u in P[1] };
> E2 := { [Index(u), t] : u in P[2] };
> E3 := { [Index(u), Index(v)] : u in P[1], v in P[2] };
>
> N := MultiDigraph< Order(G) | E1, E2, E3 >;
> N;
Multidigraph
Vertex  Neighbours
1       7 6 5 4 ;
2       5 4 7 6 ;
3       7 6 5 4 ;
4       9 ;
5       9 ;
6       9 ;
7       9 ;
8       1 3 2 ;
9       ;
> E := EdgeSet(N);
> for e in E do
>     AssignCapacity(~N, e, 1);
> end for;
First note that all the edges of N have been assigned capacity 1. Also, the edges of N are those connecting s to all the vertices in P[1], those connecting all the vertices of P[2] to t, and those connecting all the vertices of P[1] to all the vertices of P[2].

We construct a maximum flow F from s to t and we check that the capacity of the cut is indeed F.

> V := VertexSet(N);
> F := MaximumFlow(V!s, V!t);
> S := MinimumCut(V!s, V!t);
> F;
3
> S;
[ 8 ]
>
> c := 0;
> for u in S do
>     for v in OutNeighbours(u) do
>        if not v in S then
>           c +:= Capacity(u, v);
>           assert Capacity(u, v) eq Flow(u, v);
>        end if;
>     end for;
> end for;
> assert c eq F;
We now exhibit a matching in G. It will be given by the edges in E3 which are saturated, that is, whose flow is 1.
> M := [];
> for e in E3 do
>     u := V!e[1];
>     v := V!e[2];
>     if Flow(u, v) eq 1 then
>         E := Edges(u, v);
>         assert #E eq 1;
>         Append(~M, EndVertices(E[1]));
>     end if;
> end for;
> assert #M eq F;
> M;
[
  [ 1, 7 ],
  [ 3, 5 ],
  [ 2, 4 ]
]
Note that the statement
> assert #E eq 1;
makes sense as the network N has no parallel edges.

We end the example by showing that the capacity constraint and the skew symmetry is satisfied by all edges in the network.

> for e in EdgeSet(N) do
>     u := EndVertices(e)[1];
>     v := EndVertices(e)[2];
>
>     assert Flow(u, v) le Capacity(u, v);
>     assert Flow(u, v) eq  -Flow(v, u);
> end for;
Also, flow conservation holds for all vertices in V - {s, t}.
> s := V!s;
> t := V!t;
> for u in V do
>     if not u eq s and not u eq t then
>         f := 0;
>         for v in OutNeighbours(u) do
>             E := Edges(u, v);
>             f +:= Flow(E[1]);
>         end for;
>         for v in InNeighbours(u) do
>             E := Edges(v, u);
>             f -:= Flow(E[1]);
>         end for;
>
>         assert f eq 0;
>     end if;
> end for;
V2.28, 13 July 2023