Graph Algorithms

In the last chapter, we looked for shortest paths using an LP.

    Specialized algorithms are more efficient there

  • Dijkstra's algorithm
  • A*

Formally, a graph $G=(V,E)$ is a tuple consisting of a set of vertices $V$ and a set of edges $E$.

  • The vertices are arbitrary objects
    $V = \{A,B,C,D\}$
  • Edges are tuples of vertices
    $E = \{(A,B), (A,C), (C,B), (C,D)\}$
    • We initially consider directed graphs.
    • Undirected graphs work similarly.

To each edge $e \in E$ we can assign costs $\color{red} c_e$,
for example, $\color{red} c_{(C, D)} = 4$

In Python's NetworkX package, we can easily define any kind of graph.


          import networkx as nx

          g = nx.DiGraph()

          # g = nx.Graph() creates an undirected graph
        

When we create edges, the nodes are automatically generated


          g.add_edge("A", "B")
          g.add_edge("A", "C")
          g.add_edge("C", "B")
          g.add_edge("C", "D")
        

We can assign arbitrary attributes to an edge or a node (afterwards, or directly when creating the edge)


          g.edges["A", "B"]["weight"] = 1
          g.add_edge("A", "C", weight=1)

          g.nodes["A"]["demand"] = 3
        

There are different ways to access the different attributes of the graph.


          print(list(g.predecessors("C")))
          print(list(g.successors("C")))
        

We can plot our graph:


          pos = nx.spectral_layout(g)
          nx.draw(g, pos=pos, with_labels=True)
          nx.draw_networkx_edge_labels(g, pos=pos)
        

Layouting large graphs is not trivial!

We are looking for the shortest path (minimum cost $c$) from a node $v_\text{source} \in V$ to $v_\text{sink} \in V$.

Idea behind Dijkstra's Algorithm:
We explore one node at a time in the order of the costs starting from $v_\text{source}$.

Dijkstra's Algorithm:

  • Mark $v_\text{source}$ as "open".
  • Set $d_{v_\text{source}} \gets 0$
  • While $v_\text{sink}$ is not "explored":
    • Choose "open" node $v$ with the smallest $d_{v}$.
    • For all unexplored successors of $v'$ from $v$:
      • Mark $v'$ as "open"
      • Set $d_{v'} \gets \min(d_{v'}, d_v + c_{(v,v')})$
    • Mark $d$ as "explored".

https://commons.wikimedia.org/wiki/File:Dijkstras_progress_animation.gif

Dijkstra's algorithm, when properly implemented, has a worst-case complexity of
$O(|E| + |V| \log |V|)$

In both worst-case scenarios and in practice, it is significantly faster than using an LP solver.

In NetworkX


          import networkx as nx

          g = nx.DiGraph()
          g.add_edge("A", "B", weight=1)
          g.add_edge("A", "C", weight=1)
          g.add_edge("C", "B", weight=3)
          g.add_edge("C", "D", weight=4)

          nx.shortest_path(g, "A", "D", weight="weight")
        

Limitation of the Dijkstra Algorithm: Edge costs must not be negative.

The Bellman-Ford Algorithm can also handle negative edges.

The algorithm works similar to Dijkstra, but also performs updates for already explored nodes.

Runtime: $O(|E||V|)$ (worst case $O(|V|^3)$).

In NetworkX:


          import networkx as nx

          g = nx.DiGraph()
          g.add_edge("A", "B", weight=1)
          g.add_edge("A", "C", weight=1)
          g.add_edge("C", "B", weight=3)
          g.add_edge("C", "D", weight=4)
          g.add_edge("D", "B", weight=-5)

          nx.shortest_path(g, "A", "B", weight="weight", method="bellman-ford")
        

If there are negative cycles in the graph, the Bellman-Ford algorithm fails.

In this case, we usually look for the shortest simple path.

A simple path is a path that does not visit any node more than once.

The problem is NP-hard.

What if we have multiple possible start or target nodes?

Example: I want to drive to the nearest gas station, but I don't care which one.

Example: Travel from A to either C or D:

Are we stuck with the graph algorithm now?

Idea: Helper nodes with free edges from C and D:


Now we can search again from a source to a sink.

Similarly, we can connect multiple sources with an auxiliary source in an analog way.

What if we want to plan a layover.

Example: Drive from Bingen to Alzey via Mainz.

Example: We want to drive from A to F via D.

    Idea: Calculate a shortest path

  • from A to D
  • from D to F

What if we have multiple options for the layover?

Example: Drive from Bingen to Alzey, but stop at any gas station along the way.

Example: We want to travel from A to E and make a stop at C or D.

The shortest path to one of the waypoints is not necessarily part of the overall shortest path.

Idea: Build a new auxiliary graph.

  • Make a copy of the entire graph.
  • For each possible intermediate stop $v$, add an edge $(v, v')$ with cost $\color{red} c_{(v, v')} = 0$ from $v$ to its copy $v'$.
  • Find the shortest path from the original source to the copy of the sink.

We will receive one "floor" for the part before and one for the part after the layover.

For each layover, we need a copy of the original graph. The complexity grows by the number of layovers.

Example: We are looking for the shortest path, but the edge costs change over time.

During rush hour, it may be more sensible to take a less congested detour.

Idea: Time Expanded Network

Each node receives a copy every "time unit" (e.g. per minute)

Most of the time we want to have two additional helper nodes to represent multiple sources and sinks.

Here we get a copy of the graph for each time unit.

Tradeoff: Accuracy vs Runtime.

Next Extension: What if we want to know the shortest path from a source to every other node?

The Dijkstra algorithm also computes the shortest path to every node that is closer to the source than the target node.

$\Rightarrow$ We simply let the algorithm run until all nodes have been explored.

If we want to find the shortest path from every node to a sink, we can do this analogously by simply swapping the source and sink, reversing all edges, and reversing the paths at the end.

In NetworkX:


          nx.single_source_shortest_path(g, "A")

          nx.single_target_shortest_path(g, "D")
        

What if we search for the shortest path from every node to every other node?

Idea: We let Dijkstra run from each source once

Runtime: $O(|V| |E| + |V|^2 \log |V|)$

In NetworkX:


          list(nx.all_pairs_dijkstra(g))
        

Here we naturally have problems with negative costs again.

Running the Bellman-Ford algorithm from every source would have a runtime of $O(|V|^2 |E|)$
(in the worst case $O(|V|^4)$).

The Floyd-Warshall algorithm finds all shortest paths in runtime $O(|V|^3)$.

Can handle negative weights but not negative cycles.

In NetworkX:


          import networkx as nx

          g = nx.DiGraph()
          g.add_edge("A", "B", weight=1)
          g.add_edge("A", "C", weight=1)
          g.add_edge("C", "B", weight=3)
          g.add_edge("C", "D", weight=4)
          g.add_edge("D", "B", weight=-5)

          nx.floyd_warshall(g)
          # Or: list(nx.all_pairs_shortest_path(g))
        
Algorithm Runtime Use
Dijkstra $O(|E| + |V| \log |V|)$ Shortest path without negative edges
Bellman-Ford $O(|E||V|)$ Shortest path without negative cycles
Floyd-Warshall $O(|V|^3)$ All shortest paths without negative cycles

    Modeling Tricks:

  • Auxiliary nodes for multiple sources or sinks
  • Graph copy for intermediate stops
  • Graph copy for time-expanded networks

Next, we will look at algorithms to calculate flows.

Find the maximum flow from A to G.

The problem can be solved using the Edmonds-Karp algorithm (a variant of the Ford-Fulkerson method).

Basic idea: For a flow that has already been calculated, we take the graph with the remaining capacities.
We then search for any remaining capacity and extend the flow.

Original Network

Flow (5)

Residual Network

Flow (3)

Residual Network

Total Flow

Runtime: $O(|V|^2 |E|)$.

In practice, that is often enough.
However, there are also flow algorithms that run almost in $O(|E|)$.

In NetworkX:


          g = nx.DiGraph()

          g.add_edge("A", "B", capacity=5)
          g.add_edge("A", "C", capacity=2)
          g.add_edge("A", "D", capacity=4)
          g.add_edge("B", "E", capacity=3)
          g.add_edge("B", "F", capacity=7)
          g.add_edge("C", "F", capacity=8)
          g.add_edge("D", "F", capacity=3)
          g.add_edge("E", "G", capacity=9)
          g.add_edge("F", "G", capacity=5)

          nx.maximum_flow(g, "A", "G")
        

If we have multiple sources or sinks, we can solve this again with auxiliary nodes.

If an edge does not have a capacity restriction, we can set the capacity to $\infty$
(or float("inf")).

Next, we are looking for Min-Cost Flows, i.e., flows with minimal costs.

In addition to Capacities, we now also have Costs.

Various settings:

  • Maximum flow from A to G with minimum costs.
  • Each node has a Demand and the flow should have minimum costs while fulfilling all demands.

The two problems can be transformed into each other.

  • Determine the maximum flow (without costs) and set it as demand.
  • Introduce auxiliary nodes. Connect them to nodes that have demand (cost = 0, capacity = |demand|).

This problem can be efficiently solved using the Network Simplex.

A specialized variant of the Simplex algorithm that takes advantage of the graph structure.

In contrast to the general simplex, there are performance guarantees here.

Runtime $O(|V||E|\log(|V|)\log(|V|C))$
($C$ represents the maximum edge costs)

(approximately $O(|V||E|)$)

In NetworkX:


          g = nx.DiGraph()
          g.add_edge("A", "B", capacity=5, weight=4)
          g.add_edge("A", "C", capacity=2, weight=1)
          g.add_edge("A", "D", capacity=4, weight=2)
          g.add_edge("B", "E", capacity=3, weight=3)
          g.add_edge("B", "F", capacity=7, weight=1)
          g.add_edge("C", "F", capacity=8, weight=1)
          g.add_edge("D", "F", capacity=3, weight=3)
          g.add_edge("E", "G", capacity=9, weight=4)
          g.add_edge("F", "G", capacity=5, weight=5)
          g.nodes["A"]["demand"] = -8
          g.nodes["G"]["demand"] = 8

          nx.min_cost_flow(g)
        

Modeling tricks, such as time-expanded graphs, work analogously here.

Next, we want to know what the bottleneck is in a maximum flow.

What exactly is a bottleneck in a network?

A s-t Cut is a set of edges which, if removed, would disconnect the graph between vertex s and vertex t.

So, we are looking for a minimal s-t Cut (based on the capacities).

A minimal s-t cut arises automatically when calculating a maximum flow.

All edges of the cut have residual capacity 0.

There can be multiple minimal s-t cuts. We may need to select from the edges with residual capacity 0.

In NetworkX:


          g = nx.DiGraph()
          g.add_edge("A", "B", capacity=5)
          g.add_edge("A", "C", capacity=2)
          g.add_edge("A", "D", capacity=4)
          g.add_edge("B", "E", capacity=3)
          g.add_edge("B", "F", capacity=7)
          g.add_edge("C", "F", capacity=8)
          g.add_edge("D", "F", capacity=3)
          g.add_edge("E", "G", capacity=9)
          g.add_edge("F", "G", capacity=5)

          nx.minimum_cut(g, "A", "G")
        

Next, we will look at problems where we want to assign nodes pairwise (Matchings).

Example: Each node is a student and we want to find study groups of 2 people each.

Each node is a student, an edge between 2 students indicates that they can form a group.

The edges of the graph are undirected here.

We are looking for the largest possible quantity of edges, so that each node appears in at most one selected edge.

Example

Not every node necessarily finds a partner

Depending on the structure of the edges, the maximum matching can also be smaller than $\lfloor \frac{|V|}{2} \rfloor$.

Such a maximum matching in any graph can be found using Edmonds's Blossom Algorithm.

Idea: We start with a matching. If there is an augmenting path in the graph, we use it to increase the matching.

Augmenting path: A path that starts and ends with a node outside the matching, and alternates between using an edge outside and inside the matching.

It can be shown that there is always an augmenting path if the matching can still be improved.

The challenging part is to find this path in general.

Runtime: $O(|E||V|^2)$

There is an algorithm that has a runtime of $O(\sqrt{|V|} |E|)$.

In NetworkX:


          g = nx.Graph()

          g.add_edge("A", "B")
          g.add_edge("A", "C")
          g.add_edge("C", "D")
          g.add_edge("A", "D")

          nx.maximal_matching(g)
        

Here we have calculated a Maximum Cardinality Matching (maximum number of edges).

If the edges have different weights, we are more likely looking for a Maximum Weight Matching.

Example: Students indicate their preferences for who they would like to form a group with.

Maximum Weight Matching

In NetworkX:


          g = nx.Graph()

          g.add_edge("A", "B", weight=1)
          g.add_edge("A", "C", weight=2)
          g.add_edge("C", "D", weight=1)
          g.add_edge("A", "D", weight=4)
          g.add_edge("B", "D", weight=1)

          nx.max_weight_matching(g)
        

Maximum Cardinality Maximum Weight Matching

In NetworkX:


          g = nx.Graph()

          g.add_edge("A", "B", weight=1)
          g.add_edge("A", "C", weight=2)
          g.add_edge("C", "D", weight=1)
          g.add_edge("A", "D", weight=4)
          g.add_edge("B", "D", weight=1)

          nx.max_weight_matching(g, maxcardinality=True)
        

Maximum Weight (with or without Maximum Cardinality) can also be found using Edmonds's Blossom Algorithm (in $O(|E||V|^2)$ runtime).

Very often we are looking for a matching on a bipartite graph.

A graph is bipartite if its nodes can be divided into two subsets such that there are no edges between nodes of the same subset.

We have a variety of courses, each of which should be assigned exactly one room. Some courses fit better in certain rooms.

Here a Maximum Cardinality Maximum Matching makes sense.

A Bipartite Matching can also be formulated well as a flow problem:

Many matching problems are bipartite.

NetworkX has specialized algorithms for bipartite graphs.


          g = nx.Graph()

          g.add_edge("K1", "R1", weight=1)
          g.add_edge("K1", "R2", weight=3)
          g.add_edge("K1", "R3", weight=2)
          g.add_edge("K2", "R2", weight=2)
          g.add_edge("K2", "R3", weight=1)

          nx.bipartite.maximum_matching(g, ["K1", "K2"])
        

Example: We want to assign students to tutorial groups. Each group can accommodate 3 students.

How do we model with matchings where we make pairwise assignments each time?

Idea: We take a node for each possible slot in one of the tutorial groups.

Next example: We have a series of different jobs $j_1, \ldots, j_n$. Each job requires one time unit to process.

We can process jobs in parallel.

Some jobs cannot be processed simultaneously because they require the same resource.

How many time units do we need at least to process all jobs?

  • We take a node for each job.
  • If two jobs cannot be processed in parallel, we connect them with an edge.
  • We now look for a Vertex Coloring.

A Vertex Coloring is a coloring of the nodes so that no connected nodes have the same color.

Each color corresponds to a time step.

We want to use as few colors as possible.


https://en.wikipedia.org/wiki/File:Petersen_graph_3-coloring.svg

Vertex coloring is an NP-hard problem.

A heuristic approach is Greedy Coloring:

  • For each node $v$:
    • Choose the first color not assigned to any neighbor of $v$.

The order of the nodes has a significant impact on the outcome here.

We can order the nodes so that those with the most edges come first ("largest first").

Or take the node next, whose neighbors already use the most colors ("saturation largest first").


          import networkx as nx

          g = nx.fast_gnp_random_graph(500, 0.05, 3)

          coloring = nx.coloring.greedy_coloring.greedy_color(g)
          # coloring = nx.coloring.greedy_coloring.greedy_color(g, strategy="saturation_largest_first")

          max(coloring.values())
        

Modification of our example:

All jobs have fixed start and end times. A machine cannot process two jobs simultaneously.

Interval Graph


https://en.wikipedia.org/wiki/File:Interval_graph.svg

For interval graphs, Greedy Coloring always leads to an optimal solution.

For severe problems, it is often a good idea to start with a very simple heuristic first, in order to quickly achieve a functional overall system
(agile approach)

When we have a set of nodes that are fully connected among each other, we call this a clique


https://en.wikipedia.org/wiki/File:6n-graf-clique.svg

During a coloring, each node of a clique must receive a different color.

A clique of maximum size is therefore a lower bound for the Vertex Coloring.

Finding a clique of maximum size (maximum clique) is unfortunately also NP-hard.

It is easier to find a clique that cannot be enlarged further (maximum clique).

  • Can an existing clique be extended by one of its neighbors?
    • Yes: Add nodes to the clique.
    • No: The clique is at its maximum.

In NetworkX, we can enumerate all maximal cliques with find_cliques:


          import networkx as nx

          g = nx.fast_gnp_random_graph(500, 0.2, 3)

          max(len(c) for c in nx.find_cliques(g))
        

The number of maximum cliques can increase exponentially with the number of nodes.

It highly depends on how malicious the graph is.

In a clique, every node knows every other node.

In a Connected Component, the goal is to find nodes that are somehow connected.


https://en.wikipedia.org/wiki/File:Pseudoforest.svg

Connected Components can be easily found using Depth First Search in approximately linear time.

Oftentimes we look for connected components to simplify another problem.

Perhaps we can consider each component individually and break down the problem in this way.

Or we can combine the components in a subsequent problem.

Example: We are looking for dates for several meetings

  • Anna and Bob can only attend at the same time.
  • Bob and Carmen can only attend at the same time.

Anna, Bob, and Carmen can be grouped together as one entity in the appointment scheduling problem.

Example from image processing: Color all pixels of a region.

It is crucial to determine when pixels are considered connected.


https://en.wikipedia.org/wiki/File:Screenshot-Figure_1.png

We can find the connected components using connected_components:


          import networkx as nx

          g = nx.Graph()
          g.add_edge(1, 2)
          g.add_edge(2, 3)
          g.add_edge(4, 5)
          g.add_node(6)
          g.add_edge(1, 7)

          list(nx.connected_components(g))
        

In directed graphs, we distinguish between two types of connected components:

  • Weakly connected: There exists an undirected path from every node to every other node.
  • Strongly connected: There exists a directed path from every node to every other node.

          import networkx as nx

          g = nx.DiGraph()

          g.add_edge(1, 2)
          g.add_edge(2, 3)
          g.add_edge(3, 1)
          g.add_edge(2, 4)
          g.add_edge(5, 6)

          print(list(nx.weakly_connected_components(g)))
          print(list(nx.strongly_connected_components(g)))
        

When we contract all nodes of a strongly connected component, we obtain a directed acyclic graph (DAG)

With contracted_nodes, we can merge nodes together.


          import networkx as nx

          g = nx.DiGraph()
          g.add_edge(1, 2)
          g.add_edge(2, 3)
          g.add_edge(3, 1)
          g.add_edge(2, 4)
          g.add_edge(2, 5)
          g.add_edge(5, 4)

          components = list(nx.strongly_connected_components(g))
          for component in components:
              if len(component) > 1:
                  component = list(component)
                  for i in range(len(component)-1):
                      nx.contracted_nodes(g, component[-1], component[i], self_loops=False, copy=False)

          print(list(g.edges))
        

A DAG is a directed graph that has no cycles.
Every strongly connected component has size 1.


https://commons.wikimedia.org/wiki/File:Graph_Condensation.svg

DAGs are often the foundation for scheduling problems.

Example: What is the critical path in a milestone plan?


https://en.wikipedia.org/wiki/File:Pert_chart_colored.svg

Finding a longest path is generally NP-hard.

In DAGs, this can be calculated very efficiently.

In this case, the tasks correspond to the edges and the milestones correspond to the nodes.

We can have multiple edges between the same nodes.

The appropriate class in NetworkX is MultiDiGraph


          import networkx as nx

          g = nx.MultiDiGraph()
          g.add_edge(10, 30, duration=3)
          g.add_edge(10, 20, duration=4)
          g.add_edge(20, 50, duration=3)
          g.add_edge(30, 40, duration=1)
          g.add_edge(30, 50, duration=3)
          g.add_edge(40, 50, duration=3)

          nx.dag_longest_path(g)
        

An important property of DAGs is that they have a topological sorting.

A topological sorting is an order of the nodes such that all edges point from left to right.

A topological sorting is not unique.

Each topological sort provides us with an order in which we can sequentially process tasks with dependencies, for example.


          import networkx as nx

          g = nx.MultiDiGraph()
          g.add_edge(10, 30, duration=3)
          g.add_edge(10, 20, duration=4)
          g.add_edge(20, 50, duration=3)
          g.add_edge(30, 40, duration=1)
          g.add_edge(30, 50, duration=3)
          g.add_edge(40, 50, duration=3)

          list(nx.topological_sort(g))
        

Topological sorting is possible in $O(|V| + |E|)$.

We can also use it to check whether a graph is a DAG.

Topological sorting often makes other algorithms faster (e.g. longest path and shortest path)

Next example: Find a path through the graph that visits each edge exactly once.


https://en.wikipedia.org/wiki/File:HausVomNikolaus.png

This problem is called an Euler Path or Eulerian Path

When the start and end nodes are the same, it results in an Euler Cycle (Eulerian Circuit)


          import networkx as nx

          g = nx.Graph()
          g.add_edge(1,2)
          g.add_edge(1,3)
          g.add_edge(1,4)
          g.add_edge(2,3)
          g.add_edge(2,4)
          g.add_edge(3,4)
          g.add_edge(3,5)
          g.add_edge(4,5)

          list(nx.eulerian_path(g))
        

Another well-known problem is the Traveling Salesperson Problem (TSP).

We are looking for a route through the graph that visits each node and ends at the starting node.

The TSP is NP-hard.

However, there are good heuristics that work well in practice.

For example, the Christofides algorithm

The Christofides Algorithm guarantees that the solution found is at most 50% worse than the optimum.

Runtime: $O(|V|^3)$

In NetworkX: traveling_salesman_problem

There are many other algorithms in NetworkX for various problems:

  • Minimum Spanning Tree
  • Dominating Set (NP-hard)
  • Independent Set (NP-hard)
  • Planarity
  • ...