06 - Graphs: Representation and Traversal (BFS/DFS)

A network of links, a system of dependencies between tasks, a road map: they are all graphs, and in front of each you ask yourself the same two questions. How do I write it down in memory, and how do I walk it without going round in circles forever? This article answers with two representations (adjacency list and matrix) and two traversals (breadth-first, BFS, and depth-first, DFS), and uses them to solve four concrete problems: distance in edges, connected components, execution order and cycle detection. Every algorithm is implemented, run and compared with networkx; the complexity is stated next to an operation count, never next to a time. It assumes the deque queue from article 4; the trees seen in article 5 are the special case of a graph without cycles.

What a graph is and how you write it in memory

A graph is a set of vertices and a set of edges, each of which connects two vertices. Edges can be undirected (the road between two cities can be travelled both ways) or directed (one task depends on another, not the other way round). The degree of a vertex is the number of edges touching it. I call V the number of vertices and E the number of edges: all the complexities in this article are written in V and E.

The memory choice is between two structures. The adjacency list keeps, for each vertex, the list of its neighbours; the adjacency matrix keeps a V by V table with a 1 where the edge exists. The listing builds both on random graphs with V = 1,000 and increasing density, and measures two things: how many cells or entries they occupy and how many steps the question "does the edge between u and v exist?" costs on the list.

import random

def costruisci(n, archi):
    lista = [[] for _ in range(n)]
    matrice = [[0] * n for _ in range(n)]
    for u, v in archi:
        lista[u].append(v)
        lista[v].append(u)
        matrice[u][v] = matrice[v][u] = 1
    return lista, matrice

def ha_arco_lista(lista, u, v):
    passi = 0
    for w in lista[u]:  # scans u's neighbours until it finds v
        passi += 1
        if w == v:
            return True, passi
    return False, passi

rng = random.Random(5)
print("V     E        matrix cells   list items  avg steps, edge query (list)")
for n, per_vertice in ((1000, 3), (1000, 50), (1000, 250)):
    archi = set()
    while len(archi) < n * per_vertice // 2:
        u, v = rng.sample(range(n), 2)
        archi.add((min(u, v), max(u, v)))
    lista, matrice = costruisci(n, sorted(archi))
    voci = sum(len(x) for x in lista)
    passi = sum(ha_arco_lista(lista, *rng.sample(range(n), 2))[1] for _ in range(5000)) / 5000
    print(f"{n:<5} {len(archi):<8} {n * n:<14} {voci:<11} {passi:.1f}")

# Output:
# V     E        matrix cells   list items  avg steps, edge query (list)
# 1000  1500     1000000        3000        3.0
# 1000  25000    1000000        50000       48.9
# 1000  125000   1000000        250000      218.4

The matrix occupies 1,000,000 cells in all three cases, even with 1,500 real edges: its cost depends only on V squared. The list occupies 3,000, 50,000 and 250,000 entries, two per edge in an undirected graph, so it follows E. In exchange the list answers the question about the edge by scanning the neighbours: 3.0, 48.9 and 218.4 average steps over 5,000 random pairs (the answer is almost always "no", so it scans the whole list), whereas the matrix answers with a single read. The page on graphs in Algorithms, 4th edition chooses the list as the default representation; the reason can be read in the table: real graphs are almost always sparse, with E much smaller than V squared.

The practical rule that follows: use the adjacency list when you walk the graph (BFS, DFS, most of the algorithms in this series), because scanning the neighbours of a vertex costs its degree and not V. Use the matrix when the graph is dense, or when the dominant question is "is there an edge?", or when V is small, as in the all-pairs algorithms of article 7.

BFS: level-by-level traversal and distance in edges

Breadth-first search starts from a source vertex s, visits all the neighbours at distance 1, then those at distance 2, and so on, using a queue: vertices enter in order of discovery and leave in the same order. The result is that the first path by which BFS reaches a vertex is a path with the minimum number of edges. Think of a wave in a pond: it reaches the near points first and the far ones later, and never goes back.

The code is in the graph.py module, which the later examples also use. bfs records in dist who has already been seen and at what distance; the check if v not in dist is the line that stops you going round in circles in a graph with cycles: remove it and a network with a triangle never terminates. parent stores from whom each vertex was discovered, and cammino walks it backwards to rebuild the path. deque is the right queue because popleft costs O(1), whereas list.pop(0) would shift all the elements every time, as the deque documentation recalls.

from collections import deque

def bfs(adj, s):
    """Distance in number of edges from s; also returns the parents and the edges examined."""
    dist, parent = {s: 0}, {s: None}
    coda, esaminati = deque([s]), 0
    while coda:
        u = coda.popleft()
        for v in adj[u]:
            esaminati += 1
            if v not in dist:  # the "already seen" test is what prevents infinite loops
                dist[v] = dist[u] + 1
                parent[v] = u
                coda.append(v)
    return dist, parent, esaminati

def cammino(parent, v):
    out = []
    while v is not None:
        out.append(v)
        v = parent[v]
    return out[::-1]

def dfs_ricorsiva(adj, s, visti=None, ordine=None):
    visti = set() if visti is None else visti
    ordine = [] if ordine is None else ordine
    visti.add(s)
    for v in adj[s]:
        if v not in visti:
            dfs_ricorsiva(adj, v, visti, ordine)
    ordine.append(s)  # post-order: s closes after all its descendants
    return ordine

def dfs_iterativa(adj, s):
    visti, ordine = {s}, []
    pila = [(s, iter(adj[s]))]
    while pila:
        u, it = pila[-1]
        for v in it:
            if v not in visti:
                visti.add(v)
                pila.append((v, iter(adj[v])))
                break
        else:  # no new neighbour: u is finished
            ordine.append(u)
            pila.pop()
    return ordine

def componenti(adj):
    visti, quante = set(), 0
    for s in range(len(adj)):
        if s not in visti:
            quante += 1
            visti.update(bfs(adj, s)[0])
    return quante

A small example shows the output. Six stations, six two-way links; the distance is in edges and the path is rebuilt from the parents. Note the last line: the edges examined are 12, that is the sum of the degrees, because in an undirected graph every edge appears twice, once for each endpoint.

from graph import bfs, cammino

# A network of six stations; each pair is a link in both directions
nomi = ["Lecce", "Brindisi", "Taranto", "Bari", "Foggia", "Potenza"]
archi = [(0, 1), (0, 2), (1, 3), (2, 3), (3, 4), (4, 5)]
adj = [[] for _ in nomi]
for u, v in archi:
    adj[u].append(v)
    adj[v].append(u)

dist, parent, esaminati = bfs(adj, 0)
for v in range(len(nomi)):
    print(nomi[v], dist[v], "->".join(nomi[x] for x in cammino(parent, v)))
print("edges examined:", esaminati, "= sum of degrees:", sum(len(x) for x in adj))

# Output:
# Lecce 0 Lecce
# Brindisi 1 Lecce->Brindisi
# Taranto 1 Lecce->Taranto
# Bari 2 Lecce->Brindisi->Bari
# Foggia 3 Lecce->Brindisi->Bari->Foggia
# Potenza 4 Lecce->Brindisi->Bari->Foggia->Potenza
# edges examined: 12 = sum of degrees: 12

Bari is reached in 2 edges both through Brindisi and through Taranto, and BFS keeps the first path it discovers, the one through Brindisi. When several shortest paths exist, which one is returned depends on the order of the adjacency lists, not on the algorithm.

Correctness test and step count

To verify BFS and the components I use networkx (version 3.7, installed in a trial virtual environment with pip install networkx), which here is only the independent reference: the rest of the article does not use it. The listing generates 3,000 random graphs, with isolated vertices and separate components, and checks three things: the distances equal to single_source_shortest_path_length, the number of edges examined equal to the sum of the degrees of the vertices reached, and the number of components equal to number_connected_components. It then measures the edges examined on three larger graphs.

import random
import networkx as nx
from graph import bfs, dfs_ricorsiva, dfs_iterativa, componenti

rng = random.Random(1234)
for prova in range(3000):
    n = rng.randint(1, 25)
    m = rng.randint(0, n * 2)
    G = nx.Graph()
    G.add_nodes_from(range(n))
    for _ in range(m):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v)
    adj = [list(G[u]) for u in range(n)]
    s = rng.randrange(n)
    dist, parent, esaminati = bfs(adj, s)
    ref = nx.single_source_shortest_path_length(G, s)
    assert dist == ref
    assert esaminati == sum(len(adj[u]) for u in dist)
    assert set(dfs_iterativa(adj, s)) == set(dist) == set(dfs_ricorsiva(adj, s))
    assert componenti(adj) == nx.number_connected_components(G)
print("3000 random graphs: distances, reachability and components equal to networkx", nx.__version__)

print("V      E       edges examined   V+2E")
for n, e in ((1000, 3000), (10000, 30000), (100000, 300000)):
    G = nx.gnm_random_graph(n, e, seed=1)
    adj = [list(G[u]) for u in range(n)]
    dist, _, esaminati = bfs(adj, 0)
    print(f"{n:<6} {e:<7} {esaminati:<16} {n + 2 * e}  (reached {len(dist)})")

# Output:
# 3000 random graphs: distances, reachability and components equal to networkx 3.7
# V      E       edges examined   V+2E
# 1000   3000    6000             7000  (reached 1000)
# 10000  30000   59998            70000  (reached 9974)
# 100000 300000  600000           700000  (reached 99733)

The count confirms the O(V + E) complexity stated by the Sedgewick and Wayne page, which writes "BFS takes time proportional to V + E in the worst case". The edges examined are exactly 2E when the graph is connected (6,000 on 1,000 vertices and 3,000 edges), and a little fewer when some vertices are unreachable (59,998 against 60,000, because 26 vertices are isolated or in separate components). The V term is paid in other work, which I do not count here: initializing and reading the vertices. If you want the same number for another structure, the method is this: a counter inside the innermost loop, and a comparison between the counter and the formula.

DFS: go all the way down, then come back

Depth-first search follows a path as long as it can, and only when a vertex has no new neighbours left does it step back. The route is that of someone exploring a maze keeping one hand on the wall. It has a natural recursive version, where the call stack serves as memory, and an iterative one with an explicit stack. Both are in graph.py: dfs_ricorsiva and dfs_iterativa, and they return the closing order (post-order), that is the moment at which a vertex has nothing left to explore.

The recursive version has a limit that shows up immediately in Python: the depth. On a chain graph of 5,000 vertices the recursion exceeds the interpreter's limit, which by default is 1,000 and serves to avoid overflowing the C stack (sys.setrecursionlimit). Raising the limit is a patch: it moves the crash further out, or turns it into a process halt. The iterative version keeps the stack in heap memory and has no such ceiling.

from graph import dfs_ricorsiva, dfs_iterativa

n = 5000
catena = [[i + 1] if i + 1 < n else [] for i in range(n)]  # 0 -> 1 -> ... -> 4999
try:
    dfs_ricorsiva(catena, 0)
except RecursionError as e:
    print("recursive:", type(e).__name__ + ":", e)
print("iterative: visited", len(dfs_iterativa(catena, 0)), "vertices")

# Output:
# recursive: RecursionError: maximum recursion depth exceeded
# iterative: visited 5000 vertices

The error is RecursionError, an exception and not a halt, so you catch it and show it; the iterative version visits all 5,000 vertices. In dfs_iterativa the stack holds pairs (vertex, iterator over its neighbours): the for v in it loop resumes from where it stopped, and the else branch of the for fires only when the iterator is exhausted, that is when the vertex closes. As for BFS, the same page gives the cost: DFS marks the reachable vertices in time proportional to the sum of their degrees.

Execution order and cycles with DFS

A practical problem for DFS: given tasks that depend on one another, in what order do I run them? If the dependency is an edge from "task" to "what it needs first", the order sought (topological sort) exists only if there are no cycles. DFS solves the two problems together. Every vertex has three states: white (never seen), grey (open, that is on the stack) and black (closed). An edge that reaches a grey vertex closes a cycle; if no such edge is found, the closing order, reversed, puts every task before those that require it.

import random
import networkx as nx

BIANCO, GRIGIO, NERO = 0, 1, 2

def ordina(adj):
    """Topological order by DFS, or a cycle (list of vertices) if one exists."""
    colore = [BIANCO] * len(adj)
    chiusura = []
    for radice in range(len(adj)):
        if colore[radice] != BIANCO:
            continue
        colore[radice] = GRIGIO
        pila = [(radice, iter(adj[radice]))]
        cammino = [radice]
        while pila:
            u, it = pila[-1]
            for v in it:
                if colore[v] == GRIGIO:  # edge to a vertex that is still open: cycle
                    return None, cammino[cammino.index(v):] + [v]
                if colore[v] == BIANCO:
                    colore[v] = GRIGIO
                    pila.append((v, iter(adj[v])))
                    cammino.append(v)
                    break
            else:
                colore[u] = NERO
                chiusura.append(u)
                pila.pop()
                cammino.pop()
    return chiusura[::-1], None

# edge "a -> b": a depends on b, so b must run first
compiti = {"compile": ["link"], "download": ["compile"], "test": ["compile"], "release": ["test", "link"], "link": []}
nomi = list(compiti)
idx = {n: i for i, n in enumerate(nomi)}
adj = [[idx[d] for d in compiti[n]] for n in nomi]
ordine, ciclo = ordina(adj)
print("execution order:", [nomi[i] for i in reversed(ordine)])
adj[idx["link"]].append(idx["release"])  # circular dependency
ordine, ciclo = ordina(adj)
print("cycle:", [nomi[i] for i in ciclo])

rng = random.Random(77)
cicli = 0
for prova in range(3000):
    n = rng.randint(1, 15)
    G = nx.DiGraph()
    G.add_nodes_from(range(n))
    for _ in range(rng.randint(0, n * 2)):
        u, v = rng.randrange(n), rng.randrange(n)
        if u != v:
            G.add_edge(u, v)
    adj = [list(G.successors(u)) for u in range(n)]
    ordine, ciclo = ordina(adj)
    assert (ordine is not None) == nx.is_directed_acyclic_graph(G)
    if ordine is not None:
        pos = {v: i for i, v in enumerate(ordine)}
        assert sorted(pos) == list(range(n))
        assert all(pos[u] < pos[v] for u, v in G.edges)
    else:
        cicli += 1
        assert all(G.has_edge(a, b) for a, b in zip(ciclo, ciclo[1:])) and ciclo[0] == ciclo[-1]
print("3000 digraphs: acyclicity agrees with networkx;", cicli, "with a cycle, each checked edge by edge")

# Output:
# execution order: ['link', 'compile', 'download', 'test', 'release']
# cycle: ['compile', 'link', 'release', 'test', 'compile']
# 3000 digraphs: acyclicity agrees with networkx; 1128 with a cycle, each checked edge by edge

The first example, five tasks of a release, produces an order in which link comes before compile, and release is last: every dependency is respected. The second adds the circular dependency link towards release, and the function returns the cycle, compile, link, release, test, compile, instead of an impossible order. The test on 3,000 random directed graphs agrees with networkx.is_directed_acyclic_graph in every case, and for the 1,128 graphs with a cycle the returned cycle was checked edge by edge; for the others the order respects all the edges.

Why is an edge to a grey vertex really a cycle? A grey vertex is on the stack, so it is an ancestor of the current vertex in the exploration path: the edge closes the route from the ancestor to itself. Why is an edge to a black one not? That vertex is already finished, everything it reaches has already been explored without coming back. The grey/black distinction is the subtle point, and in the undirected case it simplifies, because there it is enough not to count as a cycle the edge you have just arrived from.

Where BFS is not enough: the weights counterexample

BFS minimizes the number of edges, not the cost. If the edges have weights, the path with fewer edges can be the most expensive. The listing builds a minimal one: four vertices, a direct A-D edge of 10 and an A-B-C-D path of three edges of 1.

from graph import bfs, cammino

# direct A-D costs 10; A-B-C-D costs 1+1+1
archi = [("A", "D", 10), ("A", "B", 1), ("B", "C", 1), ("C", "D", 1)]
nomi = ["A", "B", "C", "D"]
ix = {n: i for i, n in enumerate(nomi)}
adj = [[] for _ in nomi]
peso = {}
for u, v, w in archi:
    adj[ix[u]].append(ix[v])
    adj[ix[v]].append(ix[u])
    peso[(ix[u], ix[v])] = peso[(ix[v], ix[u])] = w

def costo(p):
    return sum(peso[(a, b)] for a, b in zip(p, p[1:]))

dist, parent, _ = bfs(adj, ix["A"])
p = cammino(parent, ix["D"])
print("BFS: edges", dist[ix["D"]], "path", [nomi[i] for i in p], "cost", costo(p))
alt = [ix[x] for x in "ABCD"]
print("alternative: edges", len(alt) - 1, "path", [nomi[i] for i in alt], "cost", costo(alt))

# Output:
# BFS: edges 1 path ['A', 'D'] cost 10
# alternative: edges 3 path ['A', 'B', 'C', 'D'] cost 3

BFS chooses the direct link, 1 edge and cost 10; the three-edge path costs 3. With uniform weights, or with all weights equal to 1, the two measures coincide and BFS is enough; as soon as the weights differ another algorithm is needed. Which one, and until when it works, is the subject of article 7.

What is not here, and how to choose between BFS and DFS

Strongly connected components and bridges are not covered in this article: they require DFS with discovery times and a longer analysis, which Tarjan introduced in Depth-First Search and Linear Graph Algorithms (SIAM Journal on Computing, 1972; I read its metadata, not the text). The choice between the two traversals follows the question.

• Minimum distance in edges from a source: BFS, because it explores by levels.
• Existence of a cycle, execution order, closing of subproblems: DFS, because it tells open vertices from closed ones.
• Connected components: it makes no difference, both cost V + E; BFS avoids the recursion limit without effort.

The practical difference lies in memory: BFS keeps a whole level in the queue, which in a highly branching graph can grow to a large fraction of the vertices; DFS keeps on the stack only the current path, which can be as long as V only in a chain graph, like the one in the recursion listing.

What you can do now

You can choose between list and matrix by reading the density of the graph and the dominant question, and you can write BFS and DFS without recursion that breaks on deep graphs. You can count the edges examined and compare them with V + E, you can order tasks with dependencies and return the cycle when the order does not exist, and you can build a counterexample that shows where BFS stops being enough. The next article on shortest paths starts again from that counterexample; before it, the article on trees shows the same traversals on a structure without cycles.

The same series also includes complexity and Big-O, sorting algorithms, linear search, binary search and hash tables and the closing article on dynamic programming.