Every morning, John walks to the nice little coffee shop not far from his house. He just can’t get enough of that chai latte with an extra shot and cinnamon on top. One day, immersed in his musings about the nature of reality, he wondered, between a sip and another: “If I made a different path from home every day, how long until I would have to repeat a path?” It didn’t take too long until he realized that the answer is “infintely many” (you can always walk in circles and so on). Being a reasonable man, though, John also wants to keep his journey relatively short. No more than four blocks. Sometimes five, if the weather is nice.

You can picture the city John lives in as something like this:

A city modeled as a graph

The numbered circles are either points of interest, such as John’s house or the coffee shop, or intersections between two streets. The streets themselves are represented by the blue lines. This is, of course, a graph, in which the locations or intersections are the nodes and the streets are the edges.

So how would we solve this problem computationally?

Adjacency matrices

The first thing we need to decide is how to represent this graph to be manipulated by our program. One way is by using a square matrix like this:

\[A = \begin{bmatrix} 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 1 & 1 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0\\ \end{bmatrix}\]

The matrix \(A\) was constructed by following two simple rules:

  1. If the graph has \(n\) nodes, we will need a \(n \times n\) matrix. In this case, \(n = 12\).
  2. If there is an edge between node \(i\) and node \(j\), then the value of \(a_{ij}\) (row \(i\), column \(j\) of the matrix) will be \(1\). Otherwise, the value will be \(0\).

This is called the adjacency matrix of the graph. By the way, since all streets in this graph are two-way streets (the graph is undirected), an edge from node \(i\) to node \(j\) appears twice in the matrix: in \(a_{ij}\) and \(a_{ji}\). In other words, the adjacency matrix of an undirected graph is symmetric. If you draw a line over its diagonal, the part of the matrix to the top of the diagonal is a mirror of the part on the bottom. We won’t exploit this symmetry in what follows, but it’s a handy sanity check when building the matrix by hand — and it’s worth noting that it disappears the moment we allow one-way streets, i.e. a directed graph.

Number of paths

Before we tackle the general problem posed by John, let’s consider a simpler one. Suppose we wanted to know how many ways one can go from node 5 to node 9, crossing just one edge. That’s almost like a silly question. Just look at element \(a_{5,9}\) of the matrix and you will see that the value is \(0\). There is no single edge connecting the two. That’s what the adjacency matrix is for: to tell which paths of length 1 exist in your graph.

What about length 2? Let’s think about what that means, for a moment. You want to go from node 5 to node 9 by crossing two edges. So first you need to go from node 5 to some other node (let’s call it \(k\)) and then, from node \(k\) to node 9. Of course, there may be many different choices of \(k\). Each such choice determines a path of length 2. By looking at the graph, we can see that the only choices we can make are 7 and 10. So there are two paths of length 2 in this case.

We can express this intuition more precisely using the adjacency matrix. Remember that each row shows all outbound edges from a node whereas each column shows all inbound edges to a node. Let’s then focus only on row 5 and column 9 (our source and destination, respectively):

\[r_5 = \begin{bmatrix} 0 & 1 & 0 & 1 & 0 & 0 & \mathbf{1} & 0 & 0 & \mathbf{1} & 0 & 0\\ \end{bmatrix}\] \[c_9^T = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 1 & \mathbf{1} & 1 & 0 & \mathbf{1} & 1 & 0\\ \end{bmatrix}\]

Here I’m showing column 9 horizontally (also known as its transpose) to make the point more evident. If you think of them as two arrays of numbers and look at them position by position, you will notice that there are only two positions — highlighted in bold — where the number 1 appears in both arrays: 7 and 10, matching our visual inspection of the graph.

In general, the algorithm for counting the paths of length 2 between nodes \(i\) and \(j\) is:

counter = 0
for (k = 1 to n) do
  if (a[i, k] == 1 and a[k, j] == 1) then
    counter = counter + 1;
  end
end

which is equivalent to:

counter = 0
for (k = 1 to n) do
  counter = counter + a[i, k] * a[k, j]
end

which, in mathematical notation, becomes:

\[a^{(2)}_{ij} = \sum_{k = 1}^n a_{ik}a_{kj}\]

where \(a^{(2)}_{ij}\) is the number of paths of length 2 between \(i\) and \(j\). If you are familiar with the definition of matrix multiplication, you’ll recognize this as exactly \(A^2\), shown below:

\[A^2 = \begin{bmatrix} 3 & 0 & 0 & 0 & 2 & 2 & 1 & 1 & 0 & 0 & 0 & 0\\ 0 & 2 & 1 & 2 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 3 & 2 & 0 & 0 & 1 & 0 & 2 & 0 & 1 & 0\\ 0 & 2 & 2 & 4 & 1 & 1 & 2 & 0 & 2 & 2 & 0 & 0\\ 2 & 0 & 0 & 1 & 4 & 2 & 2 & 0 & 2 & 1 & 0 & 1\\ 2 & 0 & 0 & 1 & 2 & 4 & 2 & 2 & 1 & 2 & 1 & 0\\ 1 & 1 & 1 & 2 & 2 & 2 & 5 & 1 & 2 & 2 & 1 & 1\\ 1 & 0 & 0 & 0 & 0 & 2 & 1 & 3 & 1 & 1 & 1 & 1\\ 0 & 0 & 2 & 2 & 2 & 1 & 2 & 1 & 5 & 1 & 1 & 2\\ 0 & 1 & 0 & 2 & 1 & 2 & 2 & 1 & 1 & 4 & 2 & 0\\ 0 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 2 & 3 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 1 & 2 & 0 & 0 & 2\\ \end{bmatrix}\]

But there is no reason to stop there. If we want to know all the paths of length 3, we can use the same method again. A path of length 3 is nothing more than a path of length 2 followed by a path of length 1. And we have this information already for any pair of nodes, using \(A^2\) and \(A\), respectively. So, \(A^2A = A^3\) gives us what we want. In general, if you want to know all paths of length \(m\), then \(A^m\) is your friend.

A quick but important caveat on terminology. What \(A^m\) actually counts are walks, not * *simple paths**. A walk is free to revisit nodes and edges — in our undirected graph, John could stroll from node 5 to node 7 and immediately back to node 5, and that counts as a walk of length 2. A simple path, on the other hand, never repeats a node. So some of the routes we’re counting aren’t ones John would realistically choose to walk. If we wanted to count only the simple paths, we couldn’t rely on this neat matrix trick at all: counting simple paths is a genuinely hard problem (it’s #P-complete), with no known shortcut like matrix multiplication. Walks are what make the problem tractable, so walks are what we’ll count.

The computational matrix

Notice, however, that John doesn’t need to know the number of paths from all nodes to all nodes. Only from his house to the coffee shop, which saves us a lot of computation.

Suppose his house is node 5 and the coffee shop is node 12. Instead of computing the whole multiplication of \(A\) by itself, we only need to multiply row 5 (a \(1 \times 12\) matrix) by \(A\) (a \(12 \times 12\) matrix), obtaining another \(1 \times 12\) matrix:

\[\begin{bmatrix} 2 & 0 & 0 & 1 & 4 & 2 & 2 & 0 & 2 & 1 & 0 & 1\\ \end{bmatrix}\]

In terms of graphs, this matrix contains the number of paths of length 2 from node 5 to every other node. We can then mutiply this matrix by \(A\) again to obtain the number of paths of length 3 and so on. In theory, we can repeat this infinitely many times.

This is more than a cosmetic saving. Multiplying two \(n \times n\) matrices with the schoolbook algorithm costs \(O(n^3)\) operations, so computing \(A^m\) the naive way would cost \(O(m \cdot n^3)\). By keeping only the source row — a \(1 \times n\) vector — and multiplying it by \(A\) at each step, every iteration drops to \(O(n^2)\), for a total of \(O(m \cdot n^2)\). (If we ever did need a specific high power \(A^m\) in full, there’s another trick: exponentiation by squaring computes it in \(O(\log m)\) multiplications instead of \(m\), since \(A^8 = ((A^2)^2)^2\). We don’t need it here, because we want every intermediate power up to maxLength, not just the last one.)

In Haskell:

import Data.Matrix
import Data.Vector

iteratedMult :: Vector Int -> Matrix Int -> [Vector Int]
iteratedMult row m = fmap (getRow 1) (iterate (`multStd` m) (rowVector row))

The use of the Vector type is not strictly necessary here, but makes the type of the function iteratedMult more descriptive. This function returns a lazy infinite list of all the intermediate results of this multiplication. This is where Haskell’s laziness earns its keep: iterate describes an endless sequence, but nothing is actually computed until a consumer demands it. Since we’ll later force only a finite prefix of the list (with take), the infinite construction never runs away — only as many matrix multiplications happen as there are elements we actually look at. Each intermediate result is necessary to compute the next one, but in the end all we care about are the values at the destination. So we need a function to extract those values:

pathCounts :: Int -> [Vector Int] -> [Int]
pathCounts destination = fmap (Data.Vector.! (destination - 1))

Finally, we take only as many results from this infinite list as we need and compute the sum to get our end result. Putting everything together:

numberOfPaths :: Matrix Int -> Int -> Int -> Int -> Int
numberOfPaths adjacency source destination maxLength =
  let sourceRow = getRow source adjacency
      allProducts = iteratedMult sourceRow adjacency
      paths = Prelude.take maxLength allProducts
  in Prelude.sum (pathCounts destination paths)

Using this function to solve John’s problem:

city :: Matrix Int
city = fromLists [[0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
                 ,[1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
                 ,[1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0]
                 ,[1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0]
                 ,[0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0]
                 ,[0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0]
                 ,[0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0]
                 ,[0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0]
                 ,[0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0]
                 ,[0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1]
                 ,[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1]
                 ,[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0]]

main = print $ numberOfPaths city 5 12 5
-- result: 42

There you go! In way less than 7\(\frac{1}{2}\) million years, the computer gave us the answer to the ultimate question of John, the coffee shop and everything.

A parting thought: it’s not really about counting

There’s a deeper pattern hiding in all of this, and it’s the kind of thing functional programmers tend to enjoy. Look closely at what matrix multiplication actually does: for each pair of nodes, it combines the entries along a row and column by multiplying them pairwise and then adding the results. Nothing forces those two operations to be ordinary multiplication and addition. They just need to fit together the right way — formally, they need to form a semiring.

Swap the pair \((+, \times)\) for a different one and the very same \(A^m\) machinery computes something completely different:

  • With \((+, \times)\) — the one we used — you count walks.
  • With \((\lor, \land)\) over booleans, an entry becomes “is there any walk between these nodes?”. Iterating this gives you reachability, i.e. the transitive closure of the graph.
  • With \((\min, +)\), an entry becomes the cheapest way to combine two legs of a journey. If the edges carry distances instead of 1s, \(A^m\) now computes shortest paths — this is the algebraic heart of the Floyd–Warshall algorithm.

So the code above is really one instance of a family. In Haskell you can make this precise by abstracting the matrix multiplication over a Semiring type class and letting the choice of instance decide which problem you’re solving — counting, reachability, or shortest paths — all from the same handful of lines. The graph never changed; only the algebra did.