The article below is an early draft of an article published here: https://www.toptal.com/algorithms/maximum-flow-linear-assignment-problem
Maximum Flow and the Linear Assignment Problem
Here’s a problem: Your business assigns contractors to fulfill contracts. It’s January first of the new year, and you’d like to make some assignments for the month ahead. You look through your rosters and decide which contractors are available for a one month engagement and you look through your available contracts to see which of them are for one month long tasks. I’ll also suppose that for each contractor-contract pair you know the profit you will receive if the contractor fulfills this contract. You’d like to maximize the profit of your business for this month, how will you do it?
This blog post describes an implementation of the Kuhn-Munkres algorithm (also known as the Hungarian algorithm) which solves the linear assignment problem. The Kuhn-Munkres algorithm takes time polynomial in the size of the input graph. My implementation will make use of a subroutine. The subroutine that I use will also execute in polynomial time with respect to its input (which will also be a graph). The subroutine is an implementation of the Edmonds-Karp algorithm. The Edmonds-Karp algorithm solves the maximum flow problem. The Edmonds-Karp algorithm is itself a modification of the Ford-Fulkerson method. I’ll describe the Ford-Fulkerson method first, then I’ll state the modification which turns it into the Edmonds-Karp algorithm and why that modification is important.
The maximum flow problem itself can be described informally as the problem of moving some fluid or gas through a network of pipes from a single source to a single terminal, supposing that the pressure in the network is sufficient to ensure that the fluid or gas cannot linger in any length of pipe or pipe fitting (those places where different lengths of pipe meet).
The codes in this post are in Python3.6. I used the namedtuple class in my examples.
Note: Due to some unfortunate oversights in the software I’m using to write this post underscores ( _ ) in
look like this
, but in code mode
they look like this: G_Hello!
. I switch styles pretty often so please keep in mind that if I’m talking about
and G_Hello!
I’m talking about the same thing!
Preliminaries
The ideas needed to solve these problems arise in many mathematical and engineering disciplines, often the same concepts can be known by many names and expressed in different ways (e.g. matrices or adjacency lists). Because many of the ideas are very deep, choices are made regarding how generally these concepts will be defined for any given setting. Each source makes assumptions regarding which names should be used, how the meaning behind each name will be expressed, and the level of generality in both the definitions and the discussion.
In this post I don’t want to assume any prior knowledge beyond a little introductory set theory. I want to express as little as possible and to give as few definitions as possible, making each as specific as possible (for example I don’t want to discuss infinite sets), but I still want to have what’s needed to express an graph based implementation of the Hungarian algorithm and to say something about the runtime of that implementation. I give links for sources with more information.
The implementations in this post represent the problems as directed graphs (DiGraph).
DiGraphs
A DiGraph has two attributes, setOfNodes and setOfArcs. Both these attributes are sets (unordered collections). In the code blocks on this post I’m actually using python’s frozenset but that detail isn’t particularly important.
DiGraph = namedtuple('DiGraph', ['setOfNodes','setOfArcs'])
Nodes
A Node
is composed of two attributes:
: A unique identifier.
This means that for any two Nodes
and
,
: This represent any data object.
Node = namedtuple('Node', ['uid','datum'])
Arcs
An Arc
is composed of three attributes:
: This is a Node, as defined above.
: This is a Node, as defined above.
: This represent any data object.
Arc = namedtuple('Arc', ['fromNode','toNode','datum'])
The set of Arcs in the DiGraph represents a binary relation on the Nodes in the DiGraph. Briefly, the existence of Arc
implies that a relationship exists between
and
.
In a directed graph (as opposed to an undirected graph) the existence of a relationship between
and
does not imply that a similar relationship between
and
exists,
this is because in an undirected graph the relation being expressed is not necessarily symmetric.
For example, let the Nodes in a DiGraph represent the actions necessary for building a house, and the Arcs represent the order in which these actions can occur. The relationship that the foundation must be completed first and only afterwards can the walls be affixed to the foundation can be expressed graphically like so:

In the the Arc from Node A to Node B represents the relationship allows, so if the foundation to the house has been completed, this fact allows the walls to be affixed to that foundation.
The relationship is not symmetric because the previous fact does not suggest that affixing the walls to the foundation allows the foundation to be completed.
A = Node('A','Build Foundation.')
B = Node('B','Affix walls\nto foundation.')
AB = Arc(A,B,'Allows')
DiGraphs
I’ll use Nodes and Arcs to define a DiGraph but in the algorithms below I’ll represent each DiGraph as a dictionary, because that is convenient.
Here’s a method that can convert the graph representation above into a dictionary representation similar to this one:
def digraph_to_dict(G):
G_as_dict = dict([])
for a in G.setOfArcs:
if(a.fromNode not in G.setOfNodes):
err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
pdg(G)
raise KeyError(err_msg)
if(a.toNode not in G.setOfNodes):
err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
pdg(G)
raise KeyError(err_msg)
G_as_dict[a.fromNode] = (G_as_dict[a.fromNode].union(frozenset([a]))) if (a.fromNode in G_as_dict) else frozenset([a])
for a in G.setOfArcs:
if(a.fromNode not in G.setOfNodes):
err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
if a.toNode not in G_as_dict:
G_as_dict[a.toNode] = frozenset([])
return G_as_dict
And here’s another one that converts it into a dictionary of dictionaries, another operation that will be useful:
def digraph_to_double_dict(G):
G_as_dict = dict([])
for a in G.setOfArcs:
if(a.fromNode not in G.setOfNodes):
err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
if(a.toNode not in G.setOfNodes):
err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
if(a.fromNode not in G_as_dict):
G_as_dict[a.fromNode] = dict({a.toNode : frozenset([a])})
else:
if(a.toNode not in G_as_dict[a.fromNode]):
G_as_dict[a.fromNode][a.toNode] = frozenset([a])
else:
G_as_dict[a.fromNode][a.toNode] = G_as_dict[a.fromNode][a.toNode].union(frozenset([a]))
for a in G.setOfArcs:
if(a.fromNode not in G.setOfNodes):
err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
if a.toNode not in G_as_dict:
G_as_dict[a.toNode] = dict({})
return G_as_dict
When I want to talk about a DiGraph as represented by a dictionary, I’ll use G_as_dict
to refer to it (I’ll make clear from the context whether it is a ‘double’ layer dictionary or just a single layer), otherwise I’ll use
or
.
Sometimes it’s helpful to fetch a Node from a DiGraph
by it up through it’s uid (which is unique):
def find_node_by_uid(find_uid, G):
nodes = {n for n in G.setOfNodes if n.uid == find_uid}
if(len(nodes) != 1):
err_msg = 'Node with uid {find_uid!s} is not unique.'.format(**locals())
raise KeyError(err_msg)
return nodes.pop()
In defining graphs people sometimes use the terms Node and Vertex to refer to the same concept; the same is true of the terms Arc and Edge. Two popular graph representations in python are this one which uses dictionaries and another which uses objects to represent graphs. The representation in this post is somewhere in between these two commonly used representations.
This is my DiGraph representation. There are many like it, but this one is mine.
Walks and Paths
Let
be a finite sequence (ordered collection) of Arcs in a DiGraph
such that if
is any Arc in
except for the last, and
follows
in the sequence, then there must be a Node
in
such that
.
Starting from the first Arc in
, and continuing until the last Arc in
collect (in the order encountered) all Nodes
as defined above between each two consecutive Arcs in
. Label the ordered collection of Nodes collected during this operation
.
def path_arcs_to_nodes(s_arcs):
s_nodes = list([])
arc_it = iter(s_arcs)
step_count = 0
last = None
try:
at_end = False
last = a1 = next(arc_it)
while (not at_end):
s_nodes += [a1.fromNode]
last = a2 = next(arc_it)
step_count += 1
if(a1.toNode != a2.fromNode):
err_msg = "Error at index {step_count!s} of Arc sequence.".format(**locals())
raise ValueError(err_msg)
a1 = a2
except StopIteration as e:
at_end = True
if(last is not None):
s_nodes += [last.toNode]
return s_nodes
- If any Node appears more than once in the sequence
then call
a Walk on DiGraph
.
- Otherwise, call
a Path from list(S_Nodes)[0]
to list(S_Nodes)[-1]
on DiGraph
.
Source Node
Call Node
a Source Node in DiGraph
if
is in
and
contains no Arc
such that
.
def source_nodes(G):
to_nodes = frozenset({a.toNode for a in G.setOfArcs})
sources = G.setOfNodes.difference(to_nodes)
return sources
Terminal Node
Call Node
a Terminal Node in DiGraph
if
is in
and
contains no Arc
such that
.
def terminal_nodes(G):
from_nodes = frozenset({a.fromNode for a in G.setOfArcs})
terminals = G.setOfNodes.difference(from_nodes)
return terminals
Cuts and s-t Cuts
A Cut cut
of a connected DiGraph
is a subset of Arcs from
which partitions the set of Nodes
in
.
is connected if every Node
in
has at least one Arc
in
such that either
or
, but
.
Cut = namedtuple('Cut', ['setOfCutArcs'])
The definition above refers to a subset of Arcs, but it can also define a partition of the Nodes of
. I need some functions:
For the functions predecessors_of
and successors_of
,
is a Node in set G.setOfNodes of DiGraph
, and cut
is a Cut of
:
def cut_predecessors_of(n, cut, G):
allowed_arcs = G.setOfArcs.difference(frozenset(cut.setOfCutArcs))
predecessors = frozenset({})
previous_count = len(predecessors)
reach_fringe = frozenset({n})
keep_going = True
while( keep_going ):
reachable_from = frozenset({a.fromNode for a in allowed_arcs if (a.toNode in reach_fringe)})
reach_fringe = reachable_from
predecessors = predecessors.union(reach_fringe)
current_count = len(predecessors)
keep_going = current_count - previous_count > 0
previous_count = current_count
return predecessors
def cut_successors_of(n, cut, G):
allowed_arcs = G.setOfArcs.difference(frozenset(cut.setOfCutArcs))
successors = frozenset({})
previous_count = len(successors)
reach_fringe = frozenset({n})
keep_going = True
while( keep_going ):
reachable_from = frozenset({a.toNode for a in allowed_arcs if (a.fromNode in reach_fringe)})
reach_fringe = reachable_from
successors = successors.union(reach_fringe)
current_count = len(successors)
keep_going = current_count - previous_count > 0
previous_count = current_count
return successors
Let cut
be a Cut of DiGraph
.
def get_first_part(cut, G):
firstPartFringe = frozenset({a.fromNode for a in cut.setOfCutArcs})
firstPart = firstPartFringe
for n in firstPart:
preds = cut_predecessors_of(n,cut,G)
firstPart = firstPart.union(preds)
return firstPart
def get_second_part(cut, G):
secondPartFringe = frozenset({a.toNode for a in cut.setOfCutArcs})
secondPart = secondPartFringe
for n in secondPart:
succs= cut_successors_of(n,cut,G)
secondPart = secondPart.union(succs)
return secondPart
cut
is a Cut of DiGraph
if: (get_first_part(cut, G).union(get_second_part(cut, G)) == G.setOfNodes) and (len(get_first_part(cut, G).intersect(get_second_part(cut, G))) == 0)
cut
is called an x-y Cut if (x in get_first_part(cut, G)) and (y in get_second_part(cut, G) ) and (x != y)
. When the Node
in a x-y Cut cut
is a Source Node and Node
in the x-y Cut is a Terminal Node, then this Cut is called a s-t Cut.
STCut = namedtuple('STCut', ['s','t','cut'])
Flow Networks
I can use DiGraph
to represent a flow network.
Assign each Node
, where
is in
an
that is a FlowNodeDatum
:
FlowNodeDatum = namedtuple('FlowNodeDatum', ['flowIn','flowOut'])
Assign each Arc
, where
is in
an
that is a FlowArcDatum
,
FlowArcDatum = namedtuple('FlowArcDatum', ['capacity','flow'])
and
are positive real numbers.
and
are also positive real numbers.
For every node Node
in
:
n.datum.flowIn = sum({a.datum.flow for a in G.Arcs if a.toNode == n})
n.datum.flowOut = sum({a.datum.flow for a in G.Arcs if a.fromNode == n})
DiGraph
now represents a Flow Network.
The Flow of
refers to the
for all Arcs
in
.
Feasible Flows
Let DiGraph
represent a Flow Network.
The Flow Network represented by
has Feasible Flows if:
- For every Node
in
except for Source Nodes and Terminal Nodes :
.
- For every Arc
in
:
.
Condition 1 is called a Conservation Constraint.
Condition 2 is called a Capacity Constraint.
Cut Capacity
The Cut Capacity of an s-t Cut stCut
with Source Node
and Terminal Node
of a Flow Network represented by a DiGraph
is:
def cut_capacity(stCut, G):
part_1 = get_first_part(stCut.cut,G)
part_2 = get_second_part(stCut.cut,G)
s_part = part_1 if stCut.s in part_1 else part_2
t_part = part_1 if stCut.t in part_1 else part_2
cut_capacity = sum({a.datum.capacity for a in stCut.cut.setOfCutArcs if ( (a.fromNode in s_part) and (a.toNode in t_part) )})
return cut_capacity
Minimum Capacity Cut
Let stCut = stCut(s,t,cut)
be an s-t Cut of a Flow Network represented by a DiGraph
.
stCut
is the Minimum Capacity Cut of the Flow Network represented by
if there is no other s-t Cut stCutCompetitor
in this Flow Network such that:
cut_capacity(stCut, G) < cut_capacity(stCutCompetitor, G)
Stripping the Flows Away
I would like to refer to the DiGraph that would be the result of taking a DiGraph
and stripping away all the flow data from all the Nodes in
and also all the Arcs in
.
def strip_flows(G):
new_nodes= frozenset( (Node(n.uid, FlowNodeDatum(0.0,0.0)) for n in G.setOfNodes) )
new_arcs = frozenset([])
for a in G.setOfArcs:
new_fromNode = Node(a.fromNode.uid, FlowNodeDatum(0.0,0.0))
new_toNode = Node(a.toNode.uid, FlowNodeDatum(0.0,0.0))
new_arc = Arc(new_fromNode, new_toNode, FlowArcDatum(a.datum.capacity, 0.0))
new_arcs = new_arcs.union(frozenset([new_arc]))
return DiGraph(new_nodes, new_arcs)
Maximum Flow Problem
A Flow Network represented as a DiGraph
, a Source Node
in
and a Terminal Node
in
,
can represent a Maximum Flow Problem if:
(len(list(source_nodes(G))) == 1) and (next(iter(source_nodes(G))) == s) and (len(list(terminal_nodes(G))) == 1) and (next(iter(terminal_nodes(G))) == t)
Label this representation:
MaxFlowProblemState = namedtuple('MaxFlowProblemState', ['G','sourceNodeUid','terminalNodeUid','maxFlowProblemStateUid'])
Where sourceNodeUid
=
, terminalNodeUid
=
, and maxFlowProblemStateUid
is an identifier for the problem instance.
Maximum Flow Solution
Let maxFlowProblem
represent a Maximum Flow Problem. The solution to maxFlowProblem
can be represented by a Flow Network represented as a DiGraph
.
DiGraph
is a feasible solution to the Maximum Flow Problem on input python maxFlowProblem
if:
strip_flows(maxFlowProblem.G) == strip_flows(H)
.
is a Flow Network and has Feasible Flows.
If in addition to 1 and 2:
- There can be no other Flow Network represented by DiGraph
such that strip_flows(G) == strip_flows(K)
and find_node_by_uid(t.uid,G).flowIn < find_node_by_uid(t.uid,K).flowIn
.
Then
is also an optimal solution to maxFlowProblem
.
In other words a Feasible Maximum Flow Solution can be represented by a DiGraph:
- Which is identical to DiGraph
of the corresponding Maximum Flow Problem with the exception that the n.datum.flowIn
, n.datum.flowOut
and the a.datum.flow
of any of the Nodes and Arcs may be different.
- Represents a Flow Network that has Feasible Flows.
and it can represent an Optimal Maximum Flow Solution if additionally:
- In which the
flowIn
for the Node corresponding to the Terminal Node in the Maximum Flow Problem is as large as possible (when conditions 1 and 2 are still satisfied).
If DiGraph
represents a Feasible Maximum Flow Solution : find_node_by_uid(s.uid,H).flowOut = find_node_by_uid(t.uid,H).flowIn
this follows from the Max Flow, Min Cut Theorem (discussed below). Informally, since
is assumed to have Feasible Flows this means that Flow can neither be ‘created’ (except at Source Node
) nor ‘destroyed’ (except at Terminal Node
) while crossing any (other) Node (Conservation Constraints). So, since a Maximum Flow Problem contains only a single Source Node
and a single Terminal Node
, all flow ‘created’ at
must be ‘destroyed’ at
or the Flow Network does not have Feasible Flows (the Conservation Constraint would have been violated).
Let DiGraph
represent a Feasible Maximum Flow Solution the value above is called the s-t Flow Value of
, let mfps=MaxFlowProblemState(H,maxFlowProblem.sourceNodeUid,maxFlowProblem.terminalNodeUid,maxFlowProblem.maxFlowProblemStateUid)
this means that mfps
is a successor state of maxFlowProblem
which just means that mfps
is exacly like maxFlowProblem
with the exception that the values of a.flow
for arcs
in maxFlowProblem.setOfArcs
may be different than a.flow
for arcs
in mfps.setOfArcs
.
def get_mfps_flow(mfps):
flow_from_s = find_node_by_uid(mfps.sourceNodeUid,mfps.G).datum.flowOut
flow_to_t = find_node_by_uid(mfps.terminalNodeUid,mfps.G).datum.flowIn
if( (flow_to_t - flow_from_s) > 0):
raise Exception('Infeasible s-t flow')
return flow_to_t
Here’s a visualization of a mfps
along with it’s associated maxFlowProb
each Arc
in the image has a label, these labels are
, each Node
in the image has a label, these labels are
.

s-t Cut Flow
Let mfps
represent a MaxFlowProblemState
and let stCut
represent a Cut of mfps.G
. The Cut Flow of stCut
is defined:
def get_stcut_flow(stCut,mfps):
s = find_node_by_uid(mfps.sourceNodeUid, mfps.G)
t = find_node_by_uid(mfps.terminalNodeUid, mfps.G)
part_1 = get_first_part(stCut.cut,mfps.G)
part_2 = get_second_part(stCut.cut,mfps.G)
s_part = part_1 if s in part_1 else part_2
t_part = part_1 if t in part_1 else part_2
s_t_sum = sum(a.datum.flow for a in mfps.G.setOfArcs if (a in stCut.cut.setOfCutArcs) and (a.fromNode in s_part) )
t_s_sum = sum(a.datum.flow for a in mfps.G.setOfArcs if (a in stCut.cut.setOfCutArcs) and (a.fromNode in t_part) )
cut_flow = s_t_sum - t_s_sum
return cut_flow
s-t Cut Flow is the sum of flows from the partition containing the Source Node to the partition containing the Terminal Node minus the sum of flows from the partition containing the Terminal Node to the partition containing the Source Node.
Max Flow, Min Cut
Let maxFlowProblem
represent a Maximum Flow Problem and let the solution to maxFlowProblem
be represented by a Flow Network represented as DiGraph
.
Let minStCut
be the Minimum Capacity Cut of the Flow Network represented by maxFlowProblem.G
.
Because in the Maximum Flow Problem flow originates in only a single Source Node and terminates at a single Terminal Node and because of the Capacity Constraints and the Conservation Constraints we know that the all of the flow entering maxFlowProblem.terminalNodeUid
must cross any s-t Cut, in particular it must cross minStCut
, this means:
get_flow_value(H, maxFlowProblem) <= cut_capacity(minStCut, maxFlowProblem.G)
Solving the Maximum Flow Problem
The basic idea for solving a Maximum Flow Problem maxFlowProblem
is to start with a Maximum Flow Solution represented by DiGraph
such a starting point can be H = strip_flows(maxFlowProblem.G)
. The task is then to use
and by some greedy modification of the
values of some Arcs
in
to produce another Maximum Flow Solution represented by DiGraph
such that
doesn’t can still represent a Flow Network with Feasible Flows and get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem)
. As long as this process continues the quality (get_flow_value(K, maxFlowProblem)
) of the most recently encountered Maximum Flow Solution (
) is better than any other Maximum Flow Solution that has been found. If the process reaches a point that it know that no other improvement is possible the process can terminate and it will return the optimal Maximum Flow Solution.
The description above is general and skips many proofs such as whether such a process is possible or how long it may take, I’ll give a few more details and the algorithm.
The Max Flow, Min Cut Theorem
From the book Flows in Networks by Ford and Fulkerson the statement of the Max-Flow Min-Cut Theorem (Theorem 5.1) is:
For any network the maximal flow value from
to
is equal to the minimum cut capacity of all cuts separating
and
.
Using the definitions in this post that translates to:
The solution to a maxFlowProblem
represented by a Flow Network represented as DiGraph
is optimal if:
get_flow_value(H, maxFlowProblem) = cut_capacity(minStCut, maxFlowProblem.G)
I like this proof of the theorem and wikipedia has another one.
The Max Flow, Min Cut Theorem is used to prove the correctness and completeness of Ford-Fulkerson Method
I’ll also give some a proof of the theorem in the section after Augmenting Paths.
The Ford-Fulkerson Method and the Edmonds-Karp Algorithm
CLRS Defines the Ford-Fulkerson Method like so (section 26.2):
FORD-FULKERSON-METHOD(
,
,
):
initialize flow
to 
while there exists a augmenting path
in the residual graph
augment flow
along 
Residual Graph
The Residual Graph of a Flow Network represented as the DiGraph
can be represented as a DiGraph
:
ResidualDatum = namedtuple('ResidualDatum', ['residualCapacity','action'])
def agg_n_to_u_cap(n,u,G_as_dict):
arcs_out = G_as_dict[n]
return sum([a.datum.capacity for a in arcs_out if( (a.fromNode == n) and (a.toNode == u) ) ])
def agg_n_to_u_flow(n,u,G_as_dict):
arcs_out = G_as_dict[n]
return sum([a.datum.flow for a in arcs_out if( (a.fromNode == n) and (a.toNode == u) ) ])
def get_residual_graph_of(G):
G_as_dict = digraph_to_dict(G)
residual_nodes = G.setOfNodes
residual_arcs = frozenset([])
for n in G_as_dict:
arcs_from = G_as_dict[n]
nodes_to = frozenset([find_node_by_uid(a.toNode.uid,G) for a in arcs_from])
for u in nodes_to:
n_to_u_cap_sum = agg_n_to_u_cap(n,u,G_as_dict)
n_to_u_flow_sum = agg_n_to_u_flow(n,u,G_as_dict)
if(n_to_u_cap_sum > n_to_u_flow_sum):
flow = round(n_to_u_cap_sum - n_to_u_flow_sum, TOL)
residual_arcs = residual_arcs.union( frozenset([Arc(n,u,datum=ResidualDatum(flow, 'push'))]) )
if(n_to_u_flow_sum > 0.0):
flow = round(n_to_u_flow_sum, TOL)
residual_arcs = residual_arcs.union( frozenset([Arc(u,n,datum=ResidualDatum(flow, 'pull'))]) )
return DiGraph(residual_nodes, residual_arcs)
agg_n_to_u_cap(n,u,G_as_dict)
returns the sum of
for all Arcs in the subset of
where Arc
is in the subset if
and
.
agg_n_to_u_cap(n,u,G_as_dict)
returns the sum of
for all Arcs in the subset of
where Arc
is in the subset if
and
.
Briefly, the Residual Graph
represents certain actions which can be performed on the DiGraph
.
Each pair of Nodes
in
of the Flow Network represented by DiGraph
can generate 0, 1, or 2 Arcs in the Residual Graph
of
.
- The pair
does not generate any Arcs in
if there is no Arc
in
such that
and
.
- The pair
generates the Arc
in
where
represents an Arc labeled a Push Flow Arc from
to
a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum))
if n_to_u_cap_sum > n_to_u_flow_sum
.
- The pair
generates the Arc
in
where
represents an Arc labeled a Pull Flow Arc from
to
a = Arc(n,u,datum=ResidualNode(n_to_u_cap_sum - n_to_u_flow_sum))
if n_to_u_flow_sum > 0.0
.
- Each Push Flow Arc in
represents the action of adding a total of x <= n_to_u_cap_sum - n_to_u_flow_sum
flow to Arcs in the subset of
where Arc
is in the subset if
and
.
- Each Pull Flow Arc in
represents the action of subtracting a total of x <= n_to_u_flow_sum
flow to Arcs in the subset of
where Arc
is in the subset if
and
.
Performing an individual Push or Pull action from
on the applicable Arcs in
might generate a Flow Network without Feasible Flows because the Capacity Constraints or the Conservation Constraints might be violated in the generated Flow Network.
Here’s a visualization of the Residual Graph of the previous example visualization of a Maximum Flow Solution, in the visualization each Arc
represents a.residualCapacity
.

Augmenting Path
Let maxFlowProblem
be a MaxFlowProblem, and let G_f = get_residual_graph_of(G)
be the Residual Graph of maxFlowProblem.G
.
An Augmenting Path augmentingPath
for maxFlowProblem
is any Path from find_node_by_uid(maxFlowProblem.sourceNode,G_f)
to find_node_by_uid(maxFlowProblem.terminalNode,G_f)
.
It turns out that an augmenting Path augmentingPath
can be applied to a MaxFlowSolution represented by DiGraph
generating another MaxFlowSolution represented by DiGraph
where get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem)
if H
is not optimal. Here’s how:
def augment(augmentingPath, H):
augmentingPath = list(augmentingPath)
H_as_dict = digraph_to_dict(H)
new_nodes = frozenset({})
new_arcs = frozenset({})
visited_nodes = frozenset({})
visited_arcs = frozenset({})
bottleneck_residualCapacity = min( augmentingPath, key=lambda a: a.datum.residualCapacity ).datum.residualCapacity
for x in augmentingPath:
from_node_uid = x.fromNode.uid if x.datum.action == 'push' else x.toNode.uid
to_node_uid = x.toNode.uid if x.datum.action == 'push' else x.fromNode.uid
from_node = find_node_by_uid( from_node_uid, H )
to_node = find_node_by_uid( to_node_uid, H )
load = bottleneck_residualCapacity if x.datum.action == 'push' else -bottleneck_residualCapacity
for a in H_as_dict[from_node]:
if(a.toNode == to_node):
test_sum = a.datum.flow + load
new_flow = a.datum.flow
new_from_node_flow_out = a.fromNode.datum.flowOut
new_to_node_flow_in = a.toNode.datum.flowIn
new_from_look = {n for n in new_nodes if (n.uid == a.fromNode.uid)}
new_to_look = {n for n in new_nodes if (n.uid == a.toNode.uid) }
prev_from_node = new_from_look.pop() if (len(new_from_look)>0) else a.fromNode
prev_to_node = new_to_look.pop() if (len(new_to_look)>0) else a.toNode
new_nodes = new_nodes.difference(frozenset({prev_from_node}))
new_nodes = new_nodes.difference(frozenset({prev_to_node}))
if(test_sum > a.datum.capacity):
new_flow = a.datum.capacity
new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow + a.datum.capacity
new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow + a.datum.capacity
load = test_sum - a.datum.capacity
elif(test_sum < 0.0):
new_flow = 0.0
new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow
new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow
load = test_sum
else:
new_flow = test_sum
new_from_node_flow_out = prev_from_node.datum.flowOut - a.datum.flow + new_flow
new_to_node_flow_in = prev_to_node.datum.flowIn - a.datum.flow + new_flow
load = 0.0
new_from_node_flow_out = round(new_from_node_flow_out, TOL)
new_to_node_flow_in = round(new_to_node_flow_in, TOL)
new_flow = round(new_flow, TOL)
new_from_node = Node(prev_from_node.uid, FlowNodeDatum(prev_from_node.datum.flowIn, new_from_node_flow_out))
new_to_node = Node(prev_to_node.uid, FlowNodeDatum(new_to_node_flow_in, prev_to_node.datum.flowOut))
new_arc = Arc(new_from_node, new_to_node, FlowArcDatum(a.datum.capacity, new_flow))
visited_nodes = visited_nodes.union(frozenset({a.fromNode,a.toNode}))
visited_arcs = visited_arcs.union(frozenset({a}))
new_nodes = new_nodes.union(frozenset({new_from_node, new_to_node}))
new_arcs = new_arcs.union(frozenset({new_arc}))
not_visited_nodes = H.setOfNodes.difference(visited_nodes)
not_visited_arcs = H.setOfArcs.difference(visited_arcs)
full_new_nodes = new_nodes.union(not_visited_nodes)
full_new_arcs = new_arcs.union(not_visited_arcs)
G = DiGraph(full_new_nodes, full_new_arcs)
full_new_arcs_update = frozenset([])
for a in full_new_arcs:
new_from_node = a.fromNode
new_to_node = a.toNode
new_from_node = find_node_by_uid( a.fromNode.uid, G )
new_to_node = find_node_by_uid( a.toNode.uid, G )
full_new_arcs_update = full_new_arcs_update.union( {Arc(new_from_node, new_to_node, FlowArcDatum(a.datum.capacity, a.datum.flow))} )
G = DiGraph(full_new_nodes, full_new_arcs_update)
return G
In the above TOL
is some tolerance value for rounding the flow values in the network, this is to avoid cascading imprecision if floating point calculations. So for example I used TOL = 10
to mean round to 10 significant digits.
Let K = augment(augmentingPath, H)
, then
represents a Feasible MaxFlowSolution for maxFlowProblem
. For the statement to be true, the Flow Network represented by
must have Feasible Flows (not violate the Capacity Constraint or the Conservation Constraint.
Here’s why: In the method above, each Node added to the new Flow Network represented by DiGraph
is either an exact copy of a Node from DiGraph
, or it is a Node
which has had the same number added to it’s n.datum.flowIn
as it’s n.datum.flowOut
. This means that the Conservation Constraint is satisfied in
as long as it was satisfied in
. The Conservation Constraint is satisfied because we explicitly check that any new Arc
in the network has a.datum.flow <= a.datum.capacity
thus as long as the Arcs from the set
which were copied unmodified into
do not violate the Capacity Constraint then
does not violate the Capacity Constraint.
It’s also true that get_flow_value(H, maxFlowProblem) < get_flow_value(K, maxFlowProblem)
if H
is not optimal.
Here’s why: For an Augmenting Path augmentingPath
to exist in the DiGraph representation of the Residual Graph
of a MaxFlowProblem maxFlowProblem
then the last Arc
on augmentingPath
must be a ‘push’ Arc and it must have a.toNode == maxFlowProblem.terminalNode
. An Augmenting Path is defined as one which terminates at the Terminal Node of the MaxFlowProblem for which it is an Augmenting Path. From the definition of the Residual Graph it is clear that the last Arc in an Augmenting Path on that Residual Graph must be a ‘push’ Arc because any ‘pull’ Arc
in the Augmenting Path will have b.toNode == maxFlowProblem.terminalNode
and b.fromNode != maxFlowProblem.terminalNode
from the definition of Path. Additionally, from the definition of Path it is clear that the Terminal Node is only modified once by the augment
method. Thus augment
modifies maxFlowProblem.terminalNode.flowIn
exactly once and it increases the value of maxFlowProblem.terminalNode.flowIn
because the last Arc in the augmentingPath
must be the Arc which causes the modification in maxFlowProblem.terminalNode.flowIn
during augment
. From the definition of augment
as it applies to ‘push’ Arcs the maxFlowProblem.terminalNode.flowIn
can only be increased, not decreased.
Some Proofs from Sedgewick and Wayne
The book Algorithms, fourth edition by Robert Sedgewich and Kevin Wayne has some wonderful and short proofs (pages 892-894) that will be useful. I’ll recreate them here, though I’ll use language fitting in with previous definitions. My labels for the proofs are the same as in the Sedgewick book.
Proposition E: For any DiGraph
representing a Feasible Maximum Flow Solution to a Maximum Flow Problem maxFlowProblem
, for any stCut
get_stcut_flow(stCut,H,maxFlowProblem) = get_flow_value(H, maxFlowProblem)`.
Proof: Let stCut=stCut(maxFlowProblem.sourceNode,maxFlowProblem.terminalNode,set([a for a in H.setOfArcs if a.toNode == maxFlowProblem.terminalNode]))
. Proposition E holds for stCut
directly from the definition of s-t Flow Value. Suppose that there we wish to move some Node
from the s-partition (get_first_part(stCut.cut, G)
) and into the t-partition (get_second_part(stCut.cut, G))
, to do so we need to change stCut.cut
, which could change stcut_flow = get_stcut_flow(stCut,H,maxFlowProblem)
and invalidate Proposition E. However, let’s see how the value of stcut_flow
will change as we make this change. Node
is at equilibrium meaning that the sum of flow into Node
is equal to the sum of flow out of it (this is necessary for
to represent a Feasible Solution). Notice that all flow which is part of the stcut_flow
entering Node
enters it from the s-partition (flow entering Node
from the t-partition either directly or indirectly would not have been counted in the stcut_flow
value because it is heading the wrong direction based on the definition). Additionally, all flow exiting
will eventually (direly or indirectly) flow into the Terminal Node (proved earlier). When we move Node
into the t-partition all the flow entering
from the s-partition must be added to the new stcut_flow
value, however all flow exiting
must the be subtracted from the new stcut_flow
value, the part of the flow heading directly into the t-partition is subtracted because this flow is now internal to the new t-partition and is not counted as stcut_flow
. The part of the flow from
heading into Nodes in the s-partition must also be subtracted from stcut_flow
: After
is moved into the t-partition these flows will be directed from the t-partition and into the s-partition and so must not be accounted for in the stcut_flow
, since these flows are removed the inflow into the s-partition must be reduced by the sum of these flows, and the outflow from the s-partition into the t-partition (where all flows from s-t must end up) must be reduced by the same amount. As Node
was at equilibrium prior to the process the update will have added the same value to the new stcut_flow
value as it subtracted thus leaving Proposition E true after the update. The validity of Proposition E then follows from induction on the size of the t-partition.
Here are some example Flow Networks to help visualize the not obvious cases where Proposition E holds, in the image the red areas indicate the s-partition, the blue areas represent the t-partition, and the green Arcs indicate an s-t Cut. In the second image flow between Node A and Node B increases while the flow into Terminal Node t doesn’t change.:


Corollary: No s-t Cut Flow value can exceed the capacity of any s-t Cut.
Proposition F. (Maxflow-mincut theorem): Let
be an s-t Flow the following 3 conditions are equivalent:
- There exists an st Cut whose capacity equals the value of the flow
.
is a maxflow.
- There is no augmenting path with respect to
.
Condition 1 implies condition 2 by the corollary. Condition 2 implies condition 3 because the existence of an augmenting path implies the existence of a flow with larger values, contradicting the maximality of
. Condition 3 implies condition 1: Let
be the set of all Nodes that can be reached from
with an Augmenting Path in the Residual Graph. Let
be the remaining Arcs, then
must be in
(by our assumption). The Arcs crossing from
to
then form an stCut which contains only Arcs
where either a.datum.capacity = a.datum.flow
or a.datum.flow = 0.0
. If this were otherwise then the Nodes connected by an Arc with remaining residual capacity to
would be in the set
since there would then be an Augmenting Path from
to such a Node. The flow across the stCut is equal to the stCut’s capacity (since Arcs from
to
have flow equal to capacity) and also to the value of the stFlow (by Proposition E).
This statement of the Maxflow Mincut Theorem implies the earlier statement from Flows in Networks.
Corollary. (Integrallity property) When capacities are integers, there exists an integer-valued maxflow, and the Ford-Fulkerson algorithm finds it.
Proof: Each Augmenting Path increases the flow by a positive integer (the minimum of the unused capacities in the ‘push’ Arcs and the flows in the ‘pull’ Arcs**, all of which are always positive integers.
This justifies the Ford Fulkerson Method description from CLRS. The method is to keep finding Augmenting Paths and applying augment
to the latest maxFlowSolution
coming up with better solutions, until no more Augmenting Path meaning that the latest Maximum Flow Solution is optimal.
From Ford-Fulkerson to Edmonds-Karp
The remaining questions regarding solving Maximal Flow Problems are:
- How should Augmenting Paths be constructed?
- Will the method terminate if we use real numbers and not integers?
- How long will it take to terminate (if it does)?
The Edmonds-Karp Algorithm specifies that each Augmenting Path is constructed by a Bread First Search (BFS) of the Residual Graph, it turns out that this decision of point 1 above will also force the algorithm to terminate (point 2) and allows the asymptotic time and space complexity to be determined.
First, here’s a BFS implementation:
def bfs(sourceNode, terminalNode, G_f):
G_f_as_dict = digraph_to_dict(G_f)
parent_arcs = dict([])
visited = frozenset([])
deq = deque([sourceNode])
while len(deq) > 0:
curr = deq.popleft()
if curr == terminalNode:
break
for a in G_f_as_dict.get(curr):
if (a.toNode not in visited):
visited = visited.union(frozenset([a.toNode]))
parent_arcs[a.toNode] = a
deq.extend([a.toNode])
path = list([])
curr = terminalNode
while(curr != sourceNode):
if (curr not in parent_arcs):
err_msg = 'No augmenting path from {} to {}.'.format(sourceNode.uid, terminalNode.uid)
raise StopIteration(err_msg)
path.extend([parent_arcs[curr]])
curr = parent_arcs[curr].fromNode
path.reverse()
test = deque([path[0].fromNode])
for a in path:
if(test[-1] != a.fromNode):
err_msg = 'Bad path: {}'.format(path)
raise ValueError(err_msg)
test.extend([a.toNode])
return path
I used a deque from the python collections module.
To answer question 2 above I’ll paraphrase another proof from Sedgewick and Wayne: Proposition G. The number of Augmenting Paths needed in the Edmonds-Karp algorithm with
Nodes and
Arcs is at most
. Proof: Every Augmenting Path has a bottleneck Arc- an Arc that is deleted from the Residual Graph because it corresponds either to a ‘push’ Arc that becomes filled to capacity or a ‘pull’ Arc through which the flow becomes 0. Each time an Arc becomes a bottleneck Arc the length of any Augmenting Path through it must increase by a factor or 2. This is because each Node in a Path may appear only once or not at all (from the definition of Path) since the Paths are being explored from shortest Path to longest that means that at least one more Node must be visited by the next path that goes through the particular bottleneck Node that means an additional 2 Arcs on the path before we arrive at the Node. Since the Augmenting Path is of length at most
each Arc can be on at most
Augmenting Paths, and the total number of Augmenting Paths is at most
.
The Edmonds-Karp Algorimth executes in
. If at most
Paths will be explored during the algorithm and exploring each Path with BFS is
then the most significant term of the product and hence the asymptotic complexity is
.
Let mfp
be a MaxFlowProblemState
.
def edmonds_karp(mfp):
sid, tid = mfp.sourceNodeUid, mfp.terminalNodeUid
no_more_paths = False
loop_count = 0
while(not no_more_paths):
residual_digraph = get_residual_graph_of(mfp.G)
try:
asource = find_node_by_uid(mfp.sourceNodeUid, residual_digraph)
aterm = find_node_by_uid(mfp.terminalNodeUid, residual_digraph)
apath = bfs(asource, aterm, residual_digraph)
paint_mfp_path(mfp, loop_count, apath)
G = augment(apath, mfp.G)
s = find_node_by_uid(sid, G)
t = find_node_by_uid(tid, G)
mfp = MaxFlowProblemState(G, s.uid, t.uid, mfp.maxFlowProblemStateUid)
loop_count += 1
except StopIteration as e:
no_more_paths = True
return mfp
The version above is inefficient and has worse complexity than
since it constructs a new Maximum Flow Solution and new a Residual Graph each time (rather than modifying existing DiGraphs as the algorithm advances). To get to a true
solution the algorithm must maintain both the DiGraph representing the Maximum Flow Problem State and it’s associated Residual Graph. So the algorithm must avoid iterating over Arcs and Nodes unnecessarily and update their values and associated values in the Residual Graph only as necessary.
To write a faster Edmonds Karp algorithm I re-wrote several pieces of code from the above. I hope that going through the code which generated a new DiGraph was helpful in understanding what’s going on. In the fast algorithm I use some new tricks and python data structures that I don’t want to go over in detail. I will say that a.fromNode
and a.toNode
are now treated as strings and uids to Nodes. For this code let mfps
be a MaxFlowProblemState
import uuid
def fast_get_mfps_flow(mfps):
flow_from_s = {n for n in mfps.G.setOfNodes if n.uid == mfps.sourceNodeUid}.pop().datum.flowOut
flow_to_t = {n for n in mfps.G.setOfNodes if n.uid == mfps.terminalNodeUid}.pop().datum.flowIn
if( (flow_to_t - flow_from_s) > 0.00):
raise Exception('Infeasible s-t flow')
return flow_to_t
def fast_e_k_preprocess(G):
G = strip_flows(G)
get = dict({})
get['nodes'] = dict({})
get['node_to_node_capacity'] = dict({})
get['node_to_node_flow'] = dict({})
get['arcs'] = dict({})
get['residual_arcs'] = dict({})
for a in G.setOfArcs:
if(a.fromNode not in G.setOfNodes):
err_msg = 'There is no Node {a.fromNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
if(a.toNode not in G.setOfNodes):
err_msg = 'There is no Node {a.toNode.uid!s} to match the Arc from {a.fromNode.uid!s} to {a.toNode.uid!s}'.format(**locals())
raise KeyError(err_msg)
get['nodes'][a.fromNode.uid] = a.fromNode
get['nodes'][a.toNode.uid] = a.toNode
lark = Arc(a.fromNode.uid, a.toNode.uid, FlowArcDatumWithUid(a.datum.capacity, a.datum.flow, uuid.uuid4()))
if(a.fromNode.uid not in get['arcs']):
get['arcs'][a.fromNode.uid] = dict({a.toNode.uid : dict({lark.datum.uid : lark})})
else:
if(a.toNode.uid not in get['arcs'][a.fromNode.uid]):
get['arcs'][a.fromNode.uid][a.toNode.uid] = dict({lark.datum.uid : lark})
else:
get['arcs'][a.fromNode.uid][a.toNode.uid][lark.datum.uid] = lark
for a in G.setOfArcs:
if a.toNode.uid not in get['arcs']:
get['arcs'][a.toNode.uid] = dict({})
for n in get['nodes']:
get['residual_arcs'][n] = dict()
get['node_to_node_capacity'][n] = dict()
get['node_to_node_flow'][n] = dict()
for u in get['nodes']:
n_to_u_cap_sum = sum(a.datum.capacity for a in G.setOfArcs if (a.fromNode.uid == n) and (a.toNode.uid == u) )
n_to_u_flow_sum = sum(a.datum.flow for a in G.setOfArcs if (a.fromNode.uid == n) and (a.toNode.uid == u) )
if(n_to_u_cap_sum > n_to_u_flow_sum):
flow = round(n_to_u_cap_sum - n_to_u_flow_sum, TOL)
get['residual_arcs'][n][u] = Arc(n,u,ResidualDatum(flow, 'push'))
if(n_to_u_flow_sum > 0.0):
flow = round(n_to_u_flow_sum, TOL)
get['residual_arcs'][u][n] = Arc(u,n,ResidualDatum(flow, 'pull'))
get['node_to_node_capacity'][n][u] = n_to_u_cap_sum
get['node_to_node_flow'][n][u] = n_to_u_flow_sum
return get
def fast_bfs(sid, tid, get):
parent_of = dict([])
visited = frozenset([])
deq = deque([sid])
while len(deq) > 0:
n = deq.popleft()
if n == tid:
break
for u in get['residual_arcs'][n]:
if (u not in visited):
visited = visited.union(frozenset({u}))
parent_of[u] = get['residual_arcs'][n][u]
deq.extend([u])
path = list([])
curr = tid
while(curr != sid):
if (curr not in parent_of):
err_msg = 'No augmenting path from {} to {}.'.format(sid, curr)
raise StopIteration(err_msg)
path.extend([parent_of[curr]])
curr = parent_of[curr].fromNode
path.reverse()
return path
def fast_edmonds_karp(mfps):
sid, tid = mfps.sourceNodeUid, mfps.terminalNodeUid
get = fast_e_k_preprocess(mfps.G)
no_more_paths, loop_count = False, 0
while(not no_more_paths):
try:
apath = fast_bfs(sid, tid, get)
get = fast_augment(apath, get)
loop_count += 1
except StopIteration as e:
no_more_paths = True
nodes = frozenset(get['nodes'].values())
arcs = frozenset({})
for from_node in get['arcs']:
for to_node in get['arcs'][from_node]:
for arc in get['arcs'][from_node][to_node]:
arcs |= frozenset({get['arcs'][from_node][to_node][arc]})
G = DiGraph(nodes, arcs)
mfps = MaxFlowProblemState(G, sourceNodeUid=sid, terminalNodeUid=tid, maxFlowProblemStateUid=mfps.maxFlowProblemStateUid)
return mfps
Here’s a visualization of how this algorithm solves the example Flow Network from above. The visualization shows the steps as they are reflected in the DiGraph
representing the most up-to-date Flow Network and as they are reflected in the Residual Graph of that flow network. Augmenting Paths in the Residual Graph are shown as red paths, in the DiGraph representing the problem the set of Nodes and Arcs affected by a given Augmenting Path is highlighted in green. In each case I’ll highlight the parts of the graph that will be changed (in red or green) and then show the graph after the changes (just in black).
Here’s another visualization of how this algorithm solving a different example Flow Network. Notice that this one use real numbers and contains multiple Arcs with the same fromNode
and toNode
values.
Also notice that because Arcs with a ‘pull’ ResidualDatum may be part of the Augmenting Path, the nodes affected in the DiGraph of the Flown Network may not be on a path in
! (see the picture labeled G:15).
Residual DiGraph |
 |
 |
 |
 |
 |
[mathx_R_6 |
 |
 |
 |
 |
 |
 |
 |
 |
 |
 |
 |
 |
 |
Bipartite Graphs
Suppose we have a DiGraph
,
is Bipartite if it’s possible to partition the Nodes in
into two sets (
and
) such that for any Arc
in
it cannot be true that
in
and
in
. It also cannot be true that
in
and
in
. In other words
is Bipartite if it can be partitioned into two sets of Nodes such that every Arc must connect a Node in one set to a Node in the other set.
Testing Bipartite
Suppose we have a DiGraph
, we want to test if it is Bipartite. We can do this in
by greedy coloring the graph into two colors. First, we need to generate a new DiGraph
, this graph will have will have the same set of Nodes as
but it will have more Arcs than
. Every Arc
in
will create 2 Arcs in
, the first Arc will be identical to
, the second Arc reverses the director of
( b = Arc(a.toNode,a.fromNode,a.datum)
).
Bipartition = namedtuple('Bipartition',['firstPart', 'secondPart', 'G'])
def bipartition(G):
nodes = frozenset(G.setOfNodes
arcs = frozenset(G.setOfArcs)
arcs = arcs.union( frozenset( {Arc(a.toNode,a.fromNode,a.datum) for a in G.setOfArcs} ) )
H = DiGraph(nodes, arcs)
H_as_dict = digraph_to_dict(H)
color = dict([])
some_node = next(iter(H.setOfNodes))
deq = deque([some_node])
color[some_node] = -1
while len(deq) > 0:
curr = deq.popleft()
for a in H_as_dict.get(curr):
if (a.toNode not in color):
color[a.toNode] = -1*color[curr]
deq.extend([a.toNode])
elif(color[curr] == color[a.toNode]):
print(curr,a.toNode)
raise Exception('Not Bipartite.')
firstPart = frozenset( {n for n in color if color[n] == -1 } )
secondPart = frozenset( {n for n in color if color[n] == 1 } )
if( firstPart.union(secondPart) != G.setOfNodes ):
raise Exception('Not Bipartite.')
return Bipartition(firstPart, secondPart, G)
Matchings and Maximum Matchings
Suppose we have a DiGraph
and matching
is a subset of Arcs from
matching
is a Matching if for any two Arcs
and
in matching
: len(frozenset( {a.fromNode} ).union( {a.toNode} ).union( {b.fromNode} ).union( {b.toNode} )) == 4
. In other words, no two Arcs in a Matching share a Node. Matching matching
, is a maximum matching if there is no other Matching alt_matching
in
such that len(matching) < len(alt_matching)
, in other words matching
is a Maximum Matching if it is the largest set of Arcs from
that still satisfies the definition of Matching (the addition of any Arc not already in the matching will break the Matching definition). A Maximum Matching matching
is a Perfect Matching if every for Node
in
there exists an Arc
in matching
where a.fromNode == n or a.toNode == n
.
Maximum Bipartite Matching
A Maximum Bipartite Matching is a Maximum Matching on a DiGraph
which is Bipartite. Given that
is Bipartite the problem of finding a Maximum Bipartite Matching can be transformed into a Maximum Flow Problem solvable with the Edmonds-Karp algorithm and then the Maximum Bipartite Matching can be recovered from the solution to the Maximum Flow Problem. Let bipartition
be a Bipartition of
.
To do this I need to generate a new DiGraph (
) with some new Nodes(
) and some new Arcs(
).
contains all the Nodes in
and two more Nodes
(a Source Node) and
(a Terminal Node).
will contain one Arc for each
, if an Arc
is in
and
is in bipartition.firstPart
and
is in bipartition.secondPart
then include
in
(adding a FlowArcDatum(1,0)
). If
is in bipartition.secondPart
and
is in bipartition.firstPart
then include Arc(a.toNode,a.fromNode,FlowArcDatum(1,0))
in
. The definition of a Bipartite graph ensures that no Arc connects any Nodes where both Nodes are in the same partition.
also contains an Arc from Node
to each Node in bipartition.firstPart
. Finally,
contains an Arc each Node in bipartition.secondPart
to Node
.
for all
in
.
First partition the Nodes in
the two disjoint sets (part1
and part2
) such that no Arc in
is directed from one set to the same set (this partition is possible because
is Bipartite). Next, add all Arcs in
which are directed from part1
to part2
into
. Then create a single Source Node
and a single Terminal Node
and create some more Arcs
construct a MaxFlowProblemState
def solve_mbm( bipartition ):
s = Node(uuid.uuid4(), FlowNodeDatum(0,0))
t = Node(uuid.uuid4(), FlowNodeDatum(0,0))
translate = {}
arcs = frozenset([])
for a in bipartition.G.setOfArcs:
if ( (a.fromNode in bipartition.firstPart) and (a.toNode in bipartition.secondPart) ):
fark = Arc(a.fromNode,a.toNode,FlowArcDatum(1,0))
arcs = arcs.union({fark})
translate[frozenset({a.fromNode.uid,a.toNode.uid})] = a
elif ( (a.toNode in bipartition.firstPart) and (a.fromNode in bipartition.secondPart) ):
bark = Arc(a.toNode,a.fromNode,FlowArcDatum(1,0))
arcs = arcs.union({bark})
translate[frozenset({a.fromNode.uid,a.toNode.uid})] = a
arcs1 = frozenset( {Arc(s,n,FlowArcDatum(1,0)) for n in bipartition.firstPart } )
arcs2 = frozenset( {Arc(n,t,FlowArcDatum(1,0)) for n in bipartition.secondPart } )
arcs = arcs.union(arcs1.union(arcs2))
nodes = frozenset( {Node(n.uid,FlowNodeDatum(0,0)) for n in bipartition.G.setOfNodes} ).union({s}).union({t})
G = DiGraph(nodes, arcs)
mfp = MaxFlowProblemState(G, s.uid, t.uid, 'bipartite')
result = edmonds_karp(mfp)
lookup_set = {a for a in result.G.setOfArcs if (a.datum.flow > 0) and (a.fromNode.uid != s.uid) and (a.toNode.uid != t.uid)}
matching = {translate[frozenset({a.fromNode.uid,a.toNode.uid})] for a in lookup_set}
return matching
Minimal Node Cover
A [Node Cover] in a DiGraph
is a set of Nodes (cover
) from
such that for any Arc
of
this must be true: (a.fromNode in cover) or (a.toNode in cover)
. A Minimal Node Cover is the smallest possible set of Nodes in the graph that is still a Node Cover. Kőnig’s theorem states that in a Bipartite graph the size of the Maximum Matching on that graph is equal to the size of the Minimal Node Cover, and it suggests how the Node Cover can recovered from a Maximum Matching:
Suppose we have the Bipartition bipartition
and the Maximum Matching matching
. Define a new DiGraph
,
, the Arcs in
are the union of two sets.
The first set is Arcs
in matching
, with the change that if a.fromNode in bipartition.firstPart
and a.toNode in bipartition.secondPart
then
and
are swapped in the created Arc give such Arcs a a.datum.inMatching=True
attribute to indicate that they were derived from Arcs in a Matching.
The second set is Arcs
NOT in matching
, with the change that if a.fromNode in bipartition.secondPart
and a.toNode in bipartition.firstPart
then
and
are swapped in the created Arc (give such Arcs a a.datum.inMatching=False
attribute).
Next, run a Depth First Search (DFS) starting from each Node
in bipartition.firstPart
which is neither n == a.fromNode
nor n == a.toNodes
for any Arc
in matching
. During the DFS some Nodes are visited and some are not (store this information in a n.datum.visited
field). The Minimum Node Cover is the union of the Nodes {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) }
and the Nodes {a.fromNode for a in H.setOfArcs if (a.datum.inMatching) and (not a.toNode.datum.visited)}
.
This can be shown to lead from a Maximum Matching to a Minimal Node Cover by a proof by contradiction, take some Arc
that was supposedly not covered and consider all four cases regarding whether a.fromNode
and a.toNode
belong (whether as toNode
or fromNode
) to any Arc in Matching matching
. Each case leads to a contradiction due to the order that DFS visits Nodes and the fact that matching
is a Maximum Matching.
Suppose we have a function to execute these steps an return the set of Nodes comprising the Minimal Node Cover when given the DiGraph
, and the Maximum Matching matching
:
ArcMatchingDatum = namedtuple('ArcMatchingDatum', ['inMatching' ])
NodeMatchingDatum = namedtuple('NodeMatchingDatum', ['visited'])
def dfs_helper(snodes, G):
sniter, do_stop = iter(snodes), False
visited, visited_uids = set(), set()
while(not do_stop):
try:
stack = [ next(sniter) ]
while len(stack) > 0:
curr = stack.pop()
if curr.uid not in visited_uids:
visited = visited.union( frozenset( {Node(curr.uid, NodeMatchingDatum(False))} ) )
visited_uids = visited.union(frozenset({curr.uid}))
succ = frozenset({a.toNode for a in G.setOfArcs if a.fromNode == curr})
stack.extend( succ.difference( frozenset(visited) ) )
except StopIteration as e:
stack, do_stop = [], True
return visited
def get_min_node_cover(matching, bipartition):
nodes = frozenset( { Node(n.uid, NodeMatchingDatum(False)) for n in bipartition.G.setOfNodes} )
G = DiGraph(nodes, None)
charcs = frozenset( {a for a in matching if ( (a.fromNode in bipartition.firstPart) and (a.toNode in bipartition.secondPart) )} )
arcs0 = frozenset( { Arc(find_node_by_uid(a.toNode.uid,G), find_node_by_uid(a.fromNode.uid,G), ArcMatchingDatum(True) ) for a in charcs } )
arcs1 = frozenset( { Arc(find_node_by_uid(a.fromNode.uid,G), find_node_by_uid(a.toNode.uid,G), ArcMatchingDatum(True) ) for a in matching.difference(charcs) } )
not_matching = bipartition.G.setOfArcs.difference( matching )
charcs = frozenset( {a for a in not_matching if ( (a.fromNode in bipartition.secondPart) and (a.toNode in bipartition.firstPart) )} )
arcs2 = frozenset( { Arc(find_node_by_uid(a.toNode.uid,G), find_node_by_uid(a.fromNode.uid,G), ArcMatchingDatum(False) ) for a in charcs } )
arcs3 = frozenset( { Arc(find_node_by_uid(a.fromNode.uid,G), find_node_by_uid(a.toNode.uid,G), ArcMatchingDatum(False) ) for a in not_matching.difference(charcs) } )
arcs = arcs0.union(arcs1.union(arcs2.union(arcs3)))
G = DiGraph(nodes, arcs)
bip = Bipartition({find_node_by_uid(n.uid,G) for n in bipartition.firstPart},{find_node_by_uid(n.uid,G) for n in bipartition.secondPart},G)
match_from_nodes = frozenset({find_node_by_uid(a.fromNode.uid,G) for a in matching})
match_to_nodes = frozenset({find_node_by_uid(a.toNode.uid,G) for a in matching})
snodes = bip.firstPart.difference(match_from_nodes).difference(match_to_nodes)
visited_nodes = dfs_helper(snodes, bip.G)
not_visited_nodes = bip.G.setOfNodes.difference(visited_nodes)
H = DiGraph(visited_nodes.union(not_visited_nodes), arcs)
cover1 = frozenset( {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (a.fromNode.datum.visited) ) } )
cover2 = frozenset( {a.fromNode for a in H.setOfArcs if ( (a.datum.inMatching) and (not a.toNode.datum.visited) ) } )
min_cover_nodes = cover1.union(cover2)
true_min_cover_nodes = frozenset({find_node_by_uid(n.uid, bipartition.G) for n in min_cover_nodes})
return min_cover_nodes
The Linear Assignment Problem
The linear assignment problem consists of finding a maximum weight matching in a weighted bipartite graph. Problems like the one at the very start of this post can be expressed as a Linear Assignment Problem. Given a set of workers, a set of tasks, and a function indicating the profitability of an assignment of one worker to one task, we want to maximize the sum of all assignments that we make, this is a Linear Assignment Problem. Assume that the number of tasks and workers are equal, though I will show that this assumption is easy to remove. In the implementation I represent Arc Weights with an attribute a.datum.weight
for an Arc
.
WeightArcDatum = namedtuple('WeightArcDatum', ['weight'])
Kuhn-Munkres Algorithm
The Kuhn-Munkres Algorithm solves the The Linear Assignment Problem. A good implementation can take
time, (where
is the number of Nodes in the DiGraph representing the problem). An implementation that is easier to explain takes
(for a version which regenerates DiGraphs) and
for (for a version which maintains DiGraphs). This is similar to the two different implementations of the Edmonds-Karp algorithm.
For this description I’m only working with Complete Bipartite Graphs (those where half the Nodes are in one part of the Bipartition and the other half in the second part). In the worker, task motivation this means that there are as many workers as tasks. This seems like a significant condition (what if these sets are not equal!) but it is easy to fix this issue, I talk about how to do that in the last section.
The version of the algorithm described here uses the useful concept of zero weight arcs. Unfortunately, this concept only makes sense when we are solving a Minimization (if rather than maximizing the profits of our worker-task assignments we were instead minimizing the cost of such assignments). Fortunately, it is easy to turn a Maximum Linear Assignment Problem into a Minimum Linear Assignment Problem by setting each the Arc
weights to
where M=max({a.datum.weight for a in G.setOfArcs})
. The solution to the original Maximizing Problem will be identical to the solution Minimizing Problem after the Arc weights are changed. So for the remainder assume that we make this change.
The Kuhn-Munkres Algorithm solves minimum weight matching in a weighted bipartite graph by a sequence of Maximum Matchings in unweighted Bipartite graphs. If a we find a Perfect Matching on the DiGraph representation of the Linear Assignment Problem, and if the weight of every Arc in the Matching is zero then we have found the minimum weight matching since this matching suggests that all Nodes in the DiGraph have been Matched by an Arc with the lowest possible cost (no cost can be lower than 0, based on prior definitions). No other Arcs can be added to the Matching (because all Nodes are already matched) and no Arcs should be removed from the Matching because any possible replacement Arc will have at least as great a weight value.
If we find a Maximum Matching of the subgraph of
which contains only zero weight Arcs, and it is not a Perfect Matching we don’t have a full solution (since the Matching is not Perfect). However, we can produce a new DiGraph
by changing the weights of Arcs in
in a way that new 0-weight Arcs appear and the optimal solution of
is the same as the optimal solution of
. Since we guarantee that at least one zero weight Arc is produced at each iteration we guarantee that we will arrive at a Perfect Matching in no more than |G.setOfNodes|^{2}=N^{2} such iterations.
Suppose that in Bipartition bipartition
, bipartition.firstPart
contains Nodes representing workers, and bipartition.secondPart
represents Nodes representing tasks.
The algorithm starts by generating a new DiGraph
.
. Some Arcs in
are generated from Nodes
in bipartition.firstPart
. Each such Node
generates an Arc
in
for each Arc
in
where
or
, b=Arc(a.fromNode, a.toNode, a.datum.weight - z)
where z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) ))
.
More Arcs in
are generated from Nodes
in bipartition.secondPart
. Each such Node
generates an Arc
in
for each Arc
in
where
or
, b=Arc(a.fromNode, a.toNode, ArcWeightDatum(a.datum.weight - z))
where z=min(x.datum.weight for x in G.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) ))
.
KMA: Next, form a new DiGraph
composed of only the zero weight Arcs from
, and the Nodes incident on those Arcs. Form a bipartition
on the Nodes in
, then use solve_mbm( bipartition )
to get a Maximum Matching (matching
) on
. If matching
is a Perfect Matching in
(the Arcs in matching
are incident on all Nodes in
) then the matching
is an optimal solution to the Linear Assignment Problem.
Otherwise, if matching
is not Perfect generate the Minimal Node Cover of
using node_cover = get_min_node_cover(matching, bipartition(K))
. Next, define z=min({a.datum.weight for a in H.setOfArcs if a not in node_cover})
. Define nodes = H.setOfNodes
, arcs1 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh-z)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) and (a.toNode not in node_cover)}
, arcs2 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh)) for a in H.setOfArcs if ( (a.fromNode not in node_cover) != (a.toNode not in node_cover)}
, arcs3 = {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh+z)) for a in H.setOfArcs if ( (a.fromNode in node_cover) and (a.toNode in node_cover)}
. The !=
symbol in the previous expression acts as an XOR operator. Then arcs = arcs1.union(arcs2.union(arcs3))
. Next, H=DiGraph(nodes,arcs)
. Go back to the label KMA. The algorithm continues until a Perfect Matching is produced, this Matching is also the solution to the Linear Assignment Problem.
def kuhn_munkres( bip ):
nodes = bip.G.setOfNodes
arcs = frozenset({})
original_weights = dict({})
for a in bip.G.setOfArcs:
original_weights[a.fromNode.uid,a.toNode.uid] = a.datum.weight
for n in bip.firstPart:
z = min( {x.datum.weight for x in bip.digraph.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )} )
arcs = arcs.union( frozenset({Arc(a.fromNode, a.toNode, WeightArcDatum(a.datum.weight - z))
for a in bip.digraph.setOfArcs
if ( (a.fromNode == n) or (a.toNode == n) )}) )
for n in bip.secondPart:
z = min( {x.datum.weight for x in bip.digraph.setOfArcs if ( (x.fromNode == n) or (x.toNode == n) )} )
arcs = arcs.union( frozenset({Arc(a.fromNode, a.toNode, WeightArcDatum(a.datum.weight - z))
for a in bip.digraph.setOfArcs
if ( (a.fromNode == n) or (a.toNode == n) )}) )
H = DiGraph(nodes, arcs)
assignment, value, not_done = dict({}), 0, True
while( not_done ):
zwarcs = frozenset( {a for a in H.setOfArcs if a.datum.weight == 0} )
znodes = frozenset( {n.fromNode for n in zwarcs} ).union( frozenset( {n.toNode for n in zwarcs} ) )
K = DiGraph(znodes, zwarcs)
k_bipartition = bipartition(K)
matching = solve_mbm( k_bipartition )
mnodes = frozenset({a.fromNode for a in matching}).union(frozenset({a.toNode for a in matching}))
if( len(mnodes) == len(H.setOfNodes) ):
for a in matching:
assignment[ a.fromNode.uid ] = a.toNode.uid
value += original_weights[a.fromNode.uid,a.toNode.uid]
not_done = False
else:
node_cover = get_min_node_cover(matching, bipartition(K))
z = min( frozenset( {a.datum.weight for a in H.setOfArcs if a not in node_cover} ) )
arcs1 = frozenset( {Arc(a.fromNode,a.toNode,WeightArcDatum(a.datum.weight-z))
for a in H.setOfArcs
if ( (a.fromNode not in node_cover) and (a.toNode not in node_cover))} )
arcs2 = frozenset( {Arc(a.fromNode,a.toNode,WeightArcDatum(a.datum.weight))
for a in H.setOfArcs
if ( (a.fromNode not in node_cover) != (a.toNode not in node_cover))} )
arcs3 = frozenset( {Arc(a.fromNode,a.toNode,ArcWeightDatum(a.datum.weigh+z))
for a in H.setOfArcs
if ( (a.fromNode in node_cover) and (a.toNode in node_cover))} )
nodes = H.setOfNodes
arcs = arcs1.union(arcs2.union(arcs3))
H = DiGraph(nodes,arcs)
return value, assignment
This implementation is
because it generates a new Maximum Matching matching
at each iteration, similar to the previous two implementations of Edmonds-Karp this algorithm can be modified so that it keeps track of the matching and adapts it intelligently to each iteration. When this is done the complexity becomes
. A more advanced and more recent version of this algorithm (requiring some more advanced data structures can run in
. Details of both the simpler implementation above and about the more advanced implementation can be found at this post which motivated this blog post.
None of the operations on Arc weights modify the final assignment returned by the algorithm. Here’s why: Since our input graphs are always Complete Bipartite Graphs a solution must map each Node in one partition to another Node in the second partition, via the Arc between these two Nodes. Notice that the operations performed on the Arc weights never changes the order (ordered by weight) of the Arcs incident on any particular Node. Thus when the algorithm terminates at a Perfect Complete Bipartite Matching each Node is assigned a zero weight Arc, since the relative order of the Arcs from that Node hasn’t changed during the algorithm, and since a zero weight Arc is the cheapest possible Arc and the Perfect Complete Bipartite Matching guarantees that one such Arc exists for each Node. This means that the solution generated is indeed the same as the solution from the original Linear Assignment Problem without any modification of Arc weights.
Unbalanced Assignments
It seems like the algorithm is quite limited since as described it operates only on Complete Bipartite Graphs (those where half the Nodes are in one part of the Bipartition and the other half in the second part). In the worker, task motivation this means that there are as many workers as tasks (seems quite limiting). However, there is an easy transformation that removes this restriction. Suppose that there are fewer workers than tasks, we add some dummy workers (enough to make the resulting graph a Complete Bipartite Graph). Each dummy worker has an Arc directed from the worker to each of the tasks. Each such Arc has weight 0 (placing it in a matching gives no added profit). After this change the graph is a Complete Bipartite Graph which we can solve for. Any task assigned a dummy worker is not initiated.
Suppose that there are tasks than worker, we add some dummy tasks (enough to make the resulting graph a Complete Bipartite Graph). Each dummy task has an Arc directed from each worker to the dummy task. Each such Arc has weight 0 (placing it in a matching gives no added profit). After this change the graph is a Complete Bipartite Graph which we can solve for. Any worker assigned to dummy task is not employed during the period.
Here’s another example of this type of modification.
Linear Assignment Example
Ok, finally let’s do an example with the code I’ve been using. I’m going to modify the example problem from [here][kunmonk]. We have 3 tasks: we need to Clean The Bathroom, Sweep The Floor, and Wash The Windows. The workers available to use are Alice, Bob, Charlie, and Diane. Each of the workers gives us the wage they require per task. Here are the wages per worker:
|
Clean The Bathroom |
Sweep The Floor |
Wash The Windows |
Alice |
$2 |
$3 |
$3 |
Bob |
$3 |
$2 |
$3 |
Charlie |
$3 |
$3 |
$2 |
Diane |
$9 |
$9 |
$1 |
If we want to pay the least amount of money, but still get all the tasks done, who should do what task? Start my introducing a dummy task to make the DiGraph representing the problem Bipartite.
|
Clean The Bathroom |
Sweep The Floor |
Wash The Windows |
Do Nothing |
Alice |
$2 |
$3 |
$3 |
$0 |
Bob |
$3 |
$2 |
$3 |
$0 |
Charlie |
$3 |
$3 |
$2 |
$0 |
Diane |
$9 |
$9 |
$1 |
$0 |
Supposing that the problem is encoded in a DiGraph
, then kuhn_munkres( bipartition(G) )
will solve the problem and return the assignment. It’s easy to verify that the optimal (lowest cost) assignment will cost $5. Here’s a visualization of the solution the code above generates:

Compiled with MarkdownHere, using the MarkDownHere Cheetsheet