Push–relabel maximum flow algorithm

Template:Short description In mathematical optimization, the push–relabel algorithm (alternatively, preflow–push algorithm) is an algorithm for computing maximum flows in a flow network. The name "push–relabel" comes from the two basic operations used in the algorithm. Throughout its execution, the algorithm maintains a "preflow" and gradually converts it into a maximum flow by moving flow locally between neighboring nodes using push operations under the guidance of an admissible network maintained by relabel operations. In comparison, the Ford–Fulkerson algorithm performs global augmentations that send flow following paths from the source all the way to the sink.

The push–relabel algorithm is considered one of the most efficient maximum flow algorithms. The generic algorithm has a strongly polynomial Template:Math time complexity, which is asymptotically more efficient than the Template:Math Edmonds–Karp algorithm. Specific variants of the algorithms achieve even lower time complexities. The variant based on the highest label node selection rule has Template:Math time complexity and is generally regarded as the benchmark for maximum flow algorithms. Subcubic Template:Math time complexity can be achieved using dynamic trees, although in practice it is less efficient.Template:Reference needed

The push–relabel algorithm has been extended to compute minimum cost flows. The idea of distance labels has led to a more efficient augmenting path algorithm, which in turn can be incorporated back into the push–relabel algorithm to create a variant with even higher empirical performance.

History

The concept of a preflow was originally designed by Alexander V. Karzanov and was published in 1974 in Soviet Mathematical Dokladi 15. This pre-flow algorithm also used a push operation; however, it used distances in the auxiliary network to determine where to push the flow instead of a labeling system.

The push-relabel algorithm was designed by Andrew V. Goldberg and Robert Tarjan. The algorithm was initially presented in November 1986 in STOC '86: Proceedings of the eighteenth annual ACM symposium on Theory of computing, and then officially in October 1988 as an article in the Journal of the ACM. Both papers detail a generic form of the algorithm terminating in Template:Math along with a Template:Math sequential implementation, a Template:Math implementation using dynamic trees, and parallel/distributed implementation. As explained in, Goldberg-Tarjan introduced distance labels by incorporating them into the parallel maximum flow algorithm of Yossi Shiloach and Uzi Vishkin.

Concepts

Definitions and notations

Let:

• Template:Math be a network with capacity function Template:Math,
• Template:Math a flow network, where Template:Math and Template:Math are chosen source and sink vertices respectively,
• Template:Math denote a pre-flow in Template:Mvar,
• Template:Math denote the excess function with respect to the flow Template:Mvar, defined by Template:Math,
• Template:Math denote the residual capacity function with respect to the flow Template:Mvar, defined by Template:Math,
• Template:Math being the edges where Template:Math,

and

• Template:Math denote the residual network of Template:Mvar with respect to the flow Template:Mvar.

The push–relabel algorithm uses a nonnegative integer valid labeling function which makes use of distance labels, or heights, on nodes to determine which arcs should be selected for the push operation. This labeling function is denoted by Template:Math. This function must satisfy the following conditions in order to be considered valid:

Valid labeling:

Template:Math for all Template:Math

Source condition:

Template:Math

Sink conservation:

Template:Math

In the algorithm, the label values of Template:Mvar and Template:Mvar are fixed. Template:Math is a lower bound of the unweighted distance from Template:Mvar to Template:Mvar in Template:Math  if Template:Mvar is reachable from Template:Mvar. If Template:Mvar has been disconnected from Template:Mvar, then Template:Math is a lower bound of the unweighted distance from Template:Mvar to Template:Mvar. As a result, if a valid labeling function exists, there are no Template:Math paths in Template:Math  because no such paths can be longer than Template:Math.

An arc Template:Math  is called admissible if Template:Math. The admissible network Template:Math is composed of the set of arcs Template:Math  that are admissible. The admissible network is acyclic.

For a fixed flow Template:Math, a vertex Template:Math is called active if it has positive excess with respect to Template:Math, i.e., Template:Math.

Operations

Initialization

The algorithm starts by creating a residual graph, initializing the preflow values to zero and performing a set of saturating push operations on residual arcs Template:Math exiting the source, where Template:Math. Similarly, the labels are initialized such that the label at the source is the number of nodes in the graph, Template:Math, and all other nodes are given a label of zero. Once the initialization is complete the algorithm repeatedly performs either the push or relabel operations against active nodes until no applicable operation can be performed.

Push

The push operation applies on an admissible out-arc Template:Math of an active node Template:Mvar in Template:Math. It moves Template:Math units of flow from Template:Mvar to Template:Mvar.

push(u, v):
    assert xf[u] > 0 and 𝓁[u] == 𝓁[v] + 1
    Δ = min(xf[u], c[u][v] - f[u][v])
    f[u][v] += Δ
    f[v][u] -= Δ
    xf[u] -= Δ
    xf[v] += Δ

A push operation that causes Template:Math to reach Template:Math is called a saturating push since it uses up all the available capacity of the residual arc. Otherwise, all of the excess at the node is pushed across the residual arc. This is called an unsaturating or non-saturating push.

Relabel

The relabel operation applies on an active node Template:Mvar which is neither the source nor the sink without any admissible out-arcs in Template:Math. It modifies Template:Math to be the minimum value such that an admissible out-arc is created. Note that this always increases Template:Math and never creates a steep arc, which is an arc Template:Math such that Template:Math, and Template:Math.

relabel(u):
    assert xf[u] > 0 and 𝓁[u] <= 𝓁[v] for all v such that cf[u][v] > 0
    𝓁[u] = 1 + min(𝓁[v] for all v such that cf[u][v] > 0)

Effects of push and relabel

After a push or relabel operation, Template:Math remains a valid labeling function with respect to Template:Mvar.

For a push operation on an admissible arc Template:Math, it may add an arc Template:Math to Template:Math, where Template:Math; it may also remove the arc Template:Math from Template:Math, where it effectively removes the constraint Template:Math.

To see that a relabel operation on node Template:Mvar preserves the validity of Template:Math, notice that this is trivially guaranteed by definition for the out-arcs of u in Template:Math. For the in-arcs of Template:Mvar in Template:Math, the increased Template:Math can only satisfy the constraints less tightly, not violate them.

The generic push–relabel algorithm

The generic push–relabel algorithm is used as a proof of concept only and does not contain implementation details on how to select an active node for the push and relabel operations. This generic version of the algorithm will terminate in Template:Math.

Since Template:Math, Template:Math, and there are no paths longer than Template:Math in Template:Math, in order for Template:Math to satisfy the valid labeling condition Template:Mvar must be disconnected from Template:Mvar. At initialisation, the algorithm fulfills this requirement by creating a pre-flow Template:Mvar that saturates all out-arcs of Template:Mvar, after which Template:Math is trivially valid for all Template:Math. After initialisation, the algorithm repeatedly executes an applicable push or relabel operation until no such operations apply, at which point the pre-flow has been converted into a maximum flow.

generic-push-relabel(G, c, s, t):
    create a pre-flow f that saturates all out-arcs of s
    let 𝓁[s] = |V|
    let 𝓁[v] = 0 for all v ∈ V \ {s}
    while there is an applicable push or relabel operation do
        execute the operation

Correctness

The algorithm maintains the condition that Template:Math is a valid labeling during its execution. This can be proven true by examining the effects of the push and relabel operations on the label function Template:Math. The relabel operation increases the label value by the associated minimum plus one which will always satisfy the Template:Math constraint. The push operation can send flow from Template:Mvar to Template:Mvar if Template:Math. This may add Template:Math to Template:Math and may delete Template:Math from Template:Math. The addition of Template:Math to Template:Math will not affect the valid labeling since Template:Math. The deletion of Template:Math from Template:Math removes the corresponding constraint since the valid labeling property Template:Math only applies to residual arcs in Template:Math.

If a preflow Template:Mvar and a valid labeling Template:Math for Template:Mvar exists then there is no augmenting path from Template:Mvar to Template:Mvar in the residual graph Template:Math. This can be proven by contradiction based on inequalities which arise in the labeling function when supposing that an augmenting path does exist. If the algorithm terminates, then all nodes in Template:Math are not active. This means all Template:Math have no excess flow, and with no excess the preflow Template:Mvar obeys the flow conservation constraint and can be considered a normal flow. This flow is the maximum flow according to the max-flow min-cut theorem since there is no augmenting path from Template:Mvar to Template:Mvar.

Therefore, the algorithm will return the maximum flow upon termination.

Time complexity

In order to bound the time complexity of the algorithm, we must analyze the number of push and relabel operations which occur within the main loop. The numbers of relabel, saturating push and nonsaturating push operations are analyzed separately.

In the algorithm, the relabel operation can be performed at most Template:Math times. This is because the labeling Template:Math value for any node u can never decrease, and the maximum label value is at most Template:Math for all nodes. This means the relabel operation could potentially be performed Template:Math times for all nodes Template:Math (i.e. Template:Math). This results in a bound of Template:Math for the relabel operation.

Each saturating push on an admissible arc Template:Math removes the arc from Template:Math. For the arc to be reinserted into Template:Math for another saturating push, Template:Mvar must first be relabeled, followed by a push on the arc Template:Math, then Template:Mvar must be relabeled. In the process, Template:Math increases by at least two. Therefore, there are Template:Math saturating pushes on Template:Math, and the total number of saturating pushes is at most Template:Math. This results in a time bound of Template:Math for the saturating push operations.

Bounding the number of nonsaturating pushes can be achieved via a potential argument. We use the potential function Template:Math (i.e. Template:Math is the sum of the labels of all active nodes). It is obvious that Template:Math is Template:Math initially and stays nonnegative throughout the execution of the algorithm. Both relabels and saturating pushes can increase Template:Math. However, the value of Template:Math must be equal to 0 at termination since there cannot be any remaining active nodes at the end of the algorithm's execution. This means that over the execution of the algorithm, the nonsaturating pushes must make up the difference of the relabel and saturating push operations in order for Template:Math to terminate with a value of 0. The relabel operation can increase Template:Math by at most Template:Math. A saturating push on Template:Math activates Template:Mvar if it was inactive before the push, increasing Template:Math by at most Template:Math. Hence, the total contribution of all saturating pushes operations to Template:Math is at most Template:Math. A nonsaturating push on Template:Math always deactivates Template:Mvar, but it can also activate Template:Mvar as in a saturating push. As a result, it decreases Template:Math by at least Template:Math. Since relabels and saturating pushes increase Template:Math, the total number of nonsaturating pushes must make up the difference of Template:Math. This results in a time bound of Template:Math for the nonsaturating push operations.

In sum, the algorithm executes Template:Math relabels, Template:Math saturating pushes and Template:Math nonsaturating pushes. Data structures can be designed to pick and execute an applicable operation in Template:Math time. Therefore, the time complexity of the algorithm is Template:Math.

Example

The following is a sample execution of the generic push-relabel algorithm, as defined above, on the following simple network flow graph diagram.

In the example, the Template:Nowrap and Template:Nowrap values denote the label Template:Math and excess Template:Math, respectively, of the node during the execution of the algorithm. Each residual graph in the example only contains the residual arcs with a capacity larger than zero. Each residual graph may contain multiple iterations of the Template:Nowrap loop.

Algorithm Operation(s)	Residual Graph
Initialise the residual graph by setting the preflow to values 0 and initialising the labeling.
Initial saturating push is performed across all preflow arcs out of the source, Template:Mvar.
Node Template:Mvar is relabeled in order to push its excess flow towards the sink, Template:Mvar. The excess at Template:Mvar is then pushed to Template:Mvar then Template:Mvar in two subsequent saturating pushes; which still leaves Template:Mvar with some excess.
Once again, Template:Mvar is relabeled in order to push its excess along its last remaining positive residual (i.e. push the excess back to Template:Mvar). The node Template:Mvar is then removed from the set of active nodes.
Relabel Template:Mvar and then push its excess to Template:Mvar and Template:Mvar.
Relabel Template:Mvar and then push its excess to Template:Mvar.
Relabel Template:Mvar and then push its excess to Template:Mvar.
This leaves the node Template:Mvar as the only remaining active node, but it cannot push its excess flow towards the sink. Relabel Template:Mvar and then push its excess towards the source, Template:Mvar, via the node Template:Mvar.
Push the last bit of excess at Template:Mvar back to the source, Template:Mvar. There are no remaining active nodes. The algorithm terminates and returns the maximum flow of the graph (as seen above).

The example (but with initial flow of 0) can be run here interactively.

Practical implementations

While the generic push–relabel algorithm has Template:Math time complexity, efficient implementations achieve Template:Math or lower time complexity by enforcing appropriate rules in selecting applicable push and relabel operations. The empirical performance can be further improved by heuristics.

"Current-arc" data structure and discharge operation

The "current-arc" data structure is a mechanism for visiting the in- and out-neighbors of a node in the flow network in a static circular order. If a singly linked list of neighbors is created for a node, the data structure can be as simple as a pointer into the list that steps through the list and rewinds to the head when it runs off the end.

Based on the "current-arc" data structure, the discharge operation can be defined. A discharge operation applies on an active node and repeatedly pushes flow from the node until it becomes inactive, relabeling it as necessary to create admissible arcs in the process.

discharge(u):
    while xf[u] > 0 do
        if current-arc[u] has run off the end of neighbors[u] then
            relabel(u)
            rewind current-arc[u]
        else
            let (u, v) = current-arc[u]
            if (u, v) is admissible then
                push(u, v)
            let current-arc[u] point to the next neighbor of u

Finding the next admissible edge to push on has $O(1)$ amortized complexity. The current-arc pointer only moves to the next neighbor when the edge to the current neighbor is saturated or non-admissible, and neither of these two properties can change until the active node $u$ is relabelled. Therefore, when the pointer runs off, there are no admissible unsaturated edges and we have to relabel the active node $u$, so having moved the pointer $O(V)$ times is paid for by the $O(V)$ relabel operation.

Active node selection rules

Definition of the discharge operation reduces the push–relabel algorithm to repeatedly selecting an active node to discharge. Depending on the selection rule, the algorithm exhibits different time complexities. For the sake of brevity, we ignore Template:Mvar and Template:Mvar when referring to the nodes in the following discussion.

FIFO selection rule

The FIFO push–relabel algorithm organizes the active nodes into a queue. The initial active nodes can be inserted in arbitrary order. The algorithm always removes the node at the front of the queue for discharging. Whenever an inactive node becomes active, it is appended to the back of the queue.

The algorithm has Template:Math time complexity.

Relabel-to-front selection rule

The relabel-to-front push–relabel algorithm organizes all nodes into a linked list and maintains the invariant that the list is topologically sorted with respect to the admissible network. The algorithm scans the list from front to back and performs a discharge operation on the current node if it is active. If the node is relabeled, it is moved to the front of the list, and the scan is restarted from the front.

The algorithm also has Template:Math time complexity.

Highest label selection rule

The highest-label push–relabel algorithm organizes all nodes into buckets indexed by their labels. The algorithm always selects an active node with the largest label to discharge.

The algorithm has Template:Math time complexity. If the lowest-label selection rule is used instead, the time complexity becomes Template:Math.

Implementation techniques

Although in the description of the generic push–relabel algorithm above, Template:Math is set to zero for each node u other than Template:Mvar and Template:Mvar at the beginning, it is preferable to perform a backward breadth-first search from Template:Mvar to compute exact labels.

The algorithm is typically separated into two phases. Phase one computes a maximum pre-flow by discharging only active nodes whose labels are below Template:Mvar. Phase two converts the maximum preflow into a maximum flow by returning excess flow that cannot reach Template:Mvar to Template:Mvar. It can be shown that phase two has Template:Math time complexity regardless of the order of push and relabel operations and is therefore dominated by phase one. Alternatively, it can be implemented using flow decomposition.

Heuristics are crucial to improving the empirical performance of the algorithm. Two commonly used heuristics are the gap heuristic and the global relabeling heuristic. The gap heuristic detects gaps in the labeling function. If there is a label Template:Nowrap for which there is no node Template:Mvar such that Template:Math, then any node Template:Mvar with Template:Math has been disconnected from Template:Mvar and can be relabeled to Template:Math immediately. The global relabeling heuristic periodically performs backward breadth-first search from Template:Mvar in Template:Math to compute the exact labels of the nodes. Both heuristics skip unhelpful relabel operations, which are a bottleneck of the algorithm and contribute to the ineffectiveness of dynamic trees.

Sample implementations

#include <stdlib.h>
#include <stdio.h>

#define NODES 6
#define MIN(X,Y) ((X) < (Y) ? (X) : (Y))
#define INFINITE 10000000

void push(const int * const * C, int ** F, int *excess, int u, int v) {
  int send = MIN(excess[u], C[u][v] - F[u][v]);
  F[u][v] += send;
  F[v][u] -= send;
  excess[u] -= send;
  excess[v] += send;
}

void relabel(const int * const * C, const int * const * F, int *height, int u) {
  int v;
  int min_height = INFINITE;
  for (v = 0; v < NODES; v++) {
    if (C[u][v] - F[u][v] > 0) {
      min_height = MIN(min_height, height[v]);
      height[u] = min_height + 1;
    }
  }
};

void discharge(const int * const * C, int ** F, int *excess, int *height, int *seen, int u) {
  while (excess[u] > 0) {
    if (seen[u] < NODES) {
      int v = seen[u];
      if ((C[u][v] - F[u][v] > 0) && (height[u] > height[v])) {
        push(C, F, excess, u, v);
      } else {
        seen[u] += 1;
      }
    } else {
      relabel(C, F, height, u);
      seen[u] = 0;
    }
  }
}

void moveToFront(int i, int *A) {
  int temp = A[i];
  int n;
  for (n = i; n > 0; n--) {
    A[n] = A[n-1];
  }
  A[0] = temp;
}

int pushRelabel(const int * const * C, int ** F, int source, int sink) {
  int *excess, *height, *list, *seen, i, p;

  excess = (int *) calloc(NODES, sizeof(int));
  height = (int *) calloc(NODES, sizeof(int));
  seen = (int *) calloc(NODES, sizeof(int));

  list = (int *) calloc((NODES-2), sizeof(int));

  for (i = 0, p = 0; i < NODES; i++){
    if ((i != source) && (i != sink)) {
      list[p] = i;
      p++;
    }
  }

  height[source] = NODES;
  excess[source] = INFINITE;
  for (i = 0; i < NODES; i++)
    push(C, F, excess, source, i);

  p = 0;
  while (p < NODES - 2) {
    int u = list[p];
    int old_height = height[u];
    discharge(C, F, excess, height, seen, u);
    if (height[u] > old_height) {
      moveToFront(p, list);
      p = 0;
    } else {
      p += 1;
    }
  }
  int maxflow = 0;
  for (i = 0; i < NODES; i++)
    maxflow += F[source][i];

  free(list);

  free(seen);
  free(height);
  free(excess);

  return maxflow;
}

void printMatrix(const int * const * M) {
  int i, j;
  for (i = 0; i < NODES; i++) {
    for (j = 0; j < NODES; j++)
      printf("%d\t",M[i][j]);
    printf("\n");
  }
}

int main(void) {
  int **flow, **capacities, i;
  flow = (int **) calloc(NODES, sizeof(int*));
  capacities = (int **) calloc(NODES, sizeof(int*));
  for (i = 0; i < NODES; i++) {
    flow[i] = (int *) calloc(NODES, sizeof(int));
    capacities[i] = (int *) calloc(NODES, sizeof(int));
  }

  // Sample graph
  capacities[0][1] = 2;
  capacities[0][2] = 9;
  capacities[1][2] = 1;
  capacities[1][3] = 0;
  capacities[1][4] = 0;
  capacities[2][4] = 7;
  capacities[3][5] = 7;
  capacities[4][5] = 4;

  printf("Capacity:\n");
  printMatrix(capacities);

  printf("Max Flow:\n%d\n", pushRelabel(capacities, flow, 0, 5));

  printf("Flows:\n");
  printMatrix(flow);

  return 0;
}

def relabel_to_front(C, source: int, sink: int) -> int:
    n = len(C)  # C is the capacity matrix
    F = [[0] * n for _ in range(n)]
    # residual capacity from u to v is C[u][v] - F[u][v]

    height = [0] * n  # height of node
    excess = [0] * n  # flow into node minus flow from node
    seen   = [0] * n  # neighbours seen since last relabel
    # node "queue"
    nodelist = [i for i in range(n) if i != source and i != sink]

    def push(u, v):
        send = min(excess[u], C[u][v] - F[u][v])
        F[u][v] += send
        F[v][u] -= send
        excess[u] -= send
        excess[v] += send

    def relabel(u):
        # Find smallest new height making a push possible,
        # if such a push is possible at all.
        min_height = ∞
        for v in range(n):
            if C[u][v] - F[u][v] > 0:
                min_height = min(min_height, height[v])
                height[u] = min_height + 1

    def discharge(u):
        while excess[u] > 0:
            if seen[u] < n:  # check next neighbour
                v = seen[u]
                if C[u][v] - F[u][v] > 0 and height[u] > height[v]:
                    push(u, v)
                else:
                    seen[u] += 1
            else:  # we have checked all neighbours. must relabel
                relabel(u)
                seen[u] = 0

    height[source] = n  # longest path from source to sink is less than n long
    excess[source] = ∞  # send as much flow as possible to neighbours of source
    for v in range(n):
        push(source, v)

    p = 0
    while p < len(nodelist):
        u = nodelist[p]
        old_height = height[u]
        discharge(u)
        if height[u] > old_height:
            nodelist.insert(0, nodelist.pop(p))  # move to front of list
            p = 0  # start from front of list
        else:
            p += 1

    return sum(F[source])