Finding all simple cycles in a directed graph is a classical computer science problem, for which an elegant algorithm was published by Johnson . An optimal and asymptotically faster algorithm for undirected graphs was proposed by Birmelé et al. .
The number of cycles in a complete graph grows exponentially with the number of vertices, . Even sparse graphs, such as planar graphs, can have an exponential number of cycles . A more tractable problem that arises in real applications is that of finding all simple cycles of a given length or of length bounded by , where is typically a small integer. This has applications in the analysis of social networks, communications graphs, and financial transaction graphs, etc. Despite its several applications, to the best of our knowledge, no deterministic algorithm for this problem has been proposed to date. The closest problem that has been studied is that of finding and reporting a cycle of length . Monien  gave an algorithm for any , where is the number of vertices and is the number of edges in the graph. Alon et al.  presented an algorithm that improved this bound to , reducing the running time for small values of .
In this paper, we present a relatively simple, but fast and powerful new algorithm for finding all simple cycles of length less than or equal to in sparse directed graphs. We show that for graphs with a nearly uniform vertex degree , our algorithm takes time, where is the total number of cycles discovered, , and . While our analysis for the running time addresses only a class of sparse graphs, we provide empirical and experimental evidence of the efficiency of the algorithm for general sparse graphs as well. For the class of graphs analyzed, the running time increases linearly with the size of the graph or the number of cycles discovered because for a given and a given graph, is a constant, albeit relatively large. The algorithm is not only practical to implement, but also lends itself to easy and scalable parallelization, thus extending its suitability for real applications.
The remainder of the paper is organized as follows. We start with an overview of Johnson’s algorithm  for finding all simple cycles in Section 2, since our algorithm uses a similar overall strategy to avoid futile searches. We describe our algorithm in detail in Section 3. Sections 4 and 5 contain correctness proofs and complexity analysis, respectively. Experimental results on large graphs related to real applications are presented in Section 6. Finally, Section 7 contains concluding remarks and directions for future work.
2 Johnson’s algorithm for finding all simple cycles
Johnson’s algorithm  combines depth-first search (DFS) and backtracking with an astute strategy of blocking and unblocking vertices to prevent fruitless searches and can find all simple cycles in a directed graph in time. It remains the best-known algorithm for this problem to date. At the time of its publication, it was a significant improvement over the previous fastest algorithm , whose worst-case complexity is .
In each of its steps, Johnson’s algorithm searches for cycles that start and end at a vertex , in the strongly connected component  of the subgraph induced by vertices in . Thus, step outputs all cycles in which is the vertex with the smallest index (henceforth, referred to as the ”smallest vertex”). Such cycle searches are performed for all values of from to for which the strongly connected component has more than vertex.
As the algorithm proceeds with DFS starting at vertex , it stacks and blocks the vertices encountered along the way. The search would either lead to , in which case a cycle has been discovered, or lead to a blocked vertex. Backtracking from an unsuccessful search unstacks the vertices, but leaves them blocked. It also adds the current vertex to the Blist of the most recently unstacked vertex, thus leaving a backward trail along an unsuccessful path. Blists are data structures used to remember unsuccessful portions of a path that has been explored. The blocked vertices prevent the search from following the same unsuccessful path, or a portion of it, more than once while the vertices responsible for the failure of the path are still on the stack. When backtracking from a successful path, the algorithm unblocks each unstacked vertex, and recursively unblocks the vertices in the Blist of each vertex being unblocked. The rationale is that a vertex being unstacked was the reason for blocking the vertices in the backward trails emanating from this vertex. As soon as leaves the stack after a successful search, all the vertices that it was responsible for blocking can be considered for inclusion in a cycle in a subsequent visit via a DFS path that reaches with a different stack. For a more detailed description and analysis of the algorithm, the reader is referred to the original paper .
3 Algorithm for length-constrained simple cycles
We now present an efficient algorithm for finding all elementary cycles of length less than or equal to in a directed graph where and iff there is an edge with source and destination .
|1.||function LC_CYCLES (, )|
|2.||for = 1,|
|3.||, where and|
|if and ;|
|4.||, where is the set of all vertices reachable from|
|via BFS with levels and if , , ;|
|6.||Lock() = ;|
|7.||Blist() = ;|
|9.||blen = CYCLE_SEARCH (, , , 0);|
3.1 The outer loop
Figure 1 shows the outer loop of our algorithm. Similar to Johnson’s algorithm, our algorithm also searches for cycles whose smallest vertex is in successive subgraphs for . However, the subgraph is constructed quite differently. Let be a subgraph of induced by . Then is the subgraph of induced by all the vertices in that are reachable by breath-first search (BFS) on up to levels starting with vertex . Clearly, any vertex in that is not in cannot belong to a cycle that includes and has at most edges.
|1.||integer function CYCLE_SEARCH (, , , flen)|
|2.||blen = ;|
|3.||Lock() = flen;|
|6.||if ( then|
|7.||Output (Stack ) as a new cycle of length flen;|
|8.||blen = 1;|
|9.||else if ((flen Lock() and (flen) then|
|10.||blen = MIN (blen, 1 + CYCLE_SEARCH (, , , flen));|
|13.||if (blen ) then // True only if at least one cycle found|
|14.||RELAX_LOCKS (, , blen);|
|16.||for Adjacency() if ( Blist()) then Blist() = Blist() ;|
3.2 Cycle search with a given origin
At the heart of our algorithm is the function CYCLE_SEARCH shown in Figure 2 for finding all simple cycles of length at most that originate and terminate at vertex in subgraph . When called to explore vertex , the input parameter flen is the length of the current forward path from to . The function returns the length of the shortest path from to when the current visit to results in the detection of a cycle; otherwise, it returns .
In addition to a stack, this function uses two other key data structures.
The first one is the array Lock, which stores an integer value for each vertex and helps avoid futile searches. For any vertex that has been reached at least once and is not currently on the stack, Lock() contains one plus the maximum possible length of any path from to that has the possibility of reaching back to with at most total edges in it. This value is used to permit only those paths that reach with flen Lock() to visit it and to block all other paths (Line 9). Note that Line 9 also blocks any path that already has edges and the next edge does not lead to . As a path is being explored, Lock() is set to flen when the path reaches . If this path is eventually unsuccessful, then this value of Lock indicates that a future path reaching has the potential to succeed with its current stack only if it reaches with a smaller flen value. A cycles is detected when an edge is encountered (Line 6). At this point, blen or the backward distance from is set to 1 (Line 8). Upon return from every successful call to CYCLE_SEARCH, blen at the calling vertex is updated (Line 10) so that the shortest backward distance from for the current vertex can be returned.
The second key data structure of function CYCLE_SEARCH is the Blist associated with each vertex , which stores the indices of the source vertices of all incoming edges for which no cycle containing is found with the current state of the search stack (Line 16). It can be implemented using a data structure such as B-trees , which permit fast insertion and look-up. Recall that the role of Lock() is to allow vertex to be included only in those paths from that have the potential to terminate at within or fewer total hops (edges). When the most recent visit to results in an unsuccessful search, this is achieved by setting and leaving the value of Lock() equal to flen. When is part of a successful search, then, with the help of the Blist data structures, the function RELAX_LOCKS, called at Line 14 in CYCLE_SEARCH, appropriately sets the Lock values of and the vertices that can reach it based on the length of the current shortest path from to .
|1.||function RELAX_LOCKS (, , blen)|
|2.||if (Lock() blen+1) then|
|3.||Lock() = blen+1;|
|4.||for ( Blist())|
|5.||if ( Stack) RELAX_LOCKS (, , blen+1);|
3.3 Backtracking and relaxing Lock
Figure 3 shows the recursive procedure for reassigning the Lock values when a cycle is detected. This function is called for a vertex with a corresponding integer blen, which is the length of the currently shortest known path from to . The value of Lock() is increased to blen, if possible (Line 3). An increase in Lock() is equivalent to relaxing the conditions for to be included in a future search because this is permitted for only those paths from that reach with a length strictly less than Lock(). If blen is the length of the shortest path from to , then a new path from to with length blen or less can succeed. Relaxing Lock() can also potentially relax Lock() if edge was part of a previous unsuccessful search that led to being added to the Blist of . This process of going back to relax Lock values proceeds recursively, adding 1 to blen for each backward step (Line 5). This backward traversal is gated by the condition on Line 2 and follows only those backward tracks along which the Lock values can be increased.
3.4 Performance gains through sorting and parallelism
Note that the number of vertices and edges in subgraph in Line 3 of Function LC_CYLCES in Figure 1 monotonically decreases with . Preprocessing graph to sort and renumber its vertices in the decreasing order of the sum of their incoming and outgoing edges would ensure that the density or the average degree (ratio of number of edges to the number of vertices) of also decreases monotonically with . This can potentially improve performance in two ways. First, on an average, this is likely to reduce the size of the search space for constructing , as well as the size of . Second, this has the potential to reduce the overall amount of search in CYCLE_SEARCH. Recall that for each vertex , CYCLE_SEARCH looks for cycles starting and ending in only. All other cycles that are traversed in the process are ignored. Processing denser subgraphs (presumably with more cycles) with a root vertex of greater connectivity earlier is likely to reduce the number of valid cycles that are traversed but whose reporting is delayed until these are traversed again in the right . In Section 6, we show the improvement in running time achieved by this preprocessing step.
Function LC_CYLCES can be easily and scalably parallelized, which significantly increases the algorithm’s application potential for solving real-world problem involving very large graphs. All calls to CYCLE_SEARCH on Line 9 are independent and can be performed in parallel. The amount of work in CYCLE_SEARCH generally decreases as increases, especially if the vertices of are numbered in decreasing order of their degrees. In addition to the size of , the amount of computation in CYCLE_SEARCH also depends on the number of cycles that it detects; i.e., the number of cycles in which is the smallest vertex. Therefore, an efficient parallel implementation would need to carefully batch groups of starting vertices in order to avoid load imbalance. However, other that this relatively minor consideration, parallelization of LC_CYLCES is quite straightforward.
4 Correctness proof
The following three lemmas show that function CYCLE_SEARCH in Figure 2 finds every cycle of length less than or equal to that includes vertex in exactly once.
Function CYCLE_SEARCH does not output any cycle of length greater than .
By construction, in CYCLE_SEARCH(,,,flen), flen is the length of the path from to . Line 9 ensures that flen . Therefore, if an edge is encountered on Line 6, the length of the cycle is flen, which is . ∎
Function CYCLE_SEARCH outputs every cycle of length or less in subgraph .
Assume that Lemma 2 is false, so CYCLE_SEARCH never reaches Line 7 with and stack for some valid cycle in . Without loss of generality, also assume that among the valid cycles for which this lemma fails, is the first one.
Since is a valid cycle of length ,
Since is a cycle in , must be in the adjacency list of and if CYCLE_SEARCH(,,,) reaches Line 6 with stack , the cycle is guaranteed to be detected. Therefore, for Lemma 2 to fail, the algorithm must be blocked at some edge other than in this cycle. Without loss of generality, assume that is the first blocked edge of this cycle, where and .
Since , edge can be blocked only if
on Line 9 because flen for is . This cannot happen if this is the first visit to because all Lock values were initialized to in LC_CYCLES before the call to CYCLE_SEARCH. There are only two other possibilities for Inequality 2 to be true.
The first possibility is that the last visit to , which would have happened with a stack different from the current one, did not yield a valid cycle. In this case, Lock() would have been set equal to flen for that visit. Let’s denote that flen with . By Inequality 2, it follows that
Since is a cycle, we know that a path from to of length exists. Aslo, since was visited earlier with a different stack, a cycle other than containing exists. The length of this cycle is , or . By Inequality 3, or . Therefore, by Inequality 1, . This means that this other cycle containing meets the length constraint and is a valid cycle. Since is the first cycle to be missed by the algorithm, the previous visit to , which set Lock() to a value less than or equal to , must have produced a cycle of length . This leads to a contradiction because we are considering the scenario in which the last visit to did not yield a valid cycle. This rules out the first possibility.
The second possibility is that the last visit to did yield a valid cycle, and therefore, function RELAX_LOCKS (Figure 3) was called with as the first argument. Let be the value of blen in that call to RELAX_LOCKS. By design of RELAX_LOCKS, the following must hold upon return from that call:
By virtue of being a cycle, there is path of length from to . This path must have been successfully traversed during the last visit to because is the first edge of any valid cycle to be blocked. Therefore, Line 10 in the last call to CYCLE_SEARCH(,,,flen) will ensure that
In conclusion, no edge of the cycle can be blocked in CYCLE_SEARCH if , and this cycle is guaranteed to be reported by the algorithm. This completes the proof.
Function CYCLE_SEARCH does not output any cycle , where , more than once.
Line 5 of CYCLE_SEARCH visits each element in the adjacency list of exactly once. After edge is detected and is unstacked at the end of CYCLE_SEARCH(,,,), it can never get back to with the same stack again and the path is never explored further again. ∎
5 Time complexity
We first analyze the running time of CYCLE_SEARCH for finding all cycles starting and ending at vertex in graph . In CYCLE_SEARCH, the Lock values monotonically decrease until a valid cycle is found. Additionally, a vertex cannot be revisited without decreasing its Lock value. This ensures that no vertex is visited more than times before a valid cycle is detected. After a cycle is found, the Lock values of some vertices are relaxed or increased. During the process of relaxing the Lock values, function RELAX_LOCKS does not traverse an edge unless a Lock value can be increased. In the worst case scenario, this cannot happen to any vertex at most times. Therefore, the time associated with finding a valid cycle cannot exceed , and the time complexity of CYCLE_SEARCH for finding all valid cycles in a given subgraph is , where is the number of valid cycles in which the smallest vertex is .
Recall that consists of and all vertices greater than that are reachable from within steps of BFS. If is a sparse graph with a uniform degree , then is for all . Therefore, the complexity of the overall algorithm for finding all cycles of length or less in a graph with a uniform degree is = , where and is the total number of valid cycles in .
In our analysis, we use as the upper bound on the size of any in the call to CYCLE_SEARCH. For sufficiently large graphs, the size of would not increase with as long as stays constant. consists of all vertices with a maximum distance of from , and the size of this feasible neighborhood of depends only the structure of (unless is a scale-free network, which is rare ). The simplifying assumption of uniform vertex degrees was made to arrive at a closed-form expression for an upper bound on the size of . For sufficiently large , the size of is likely to be independent of even without this assumption. As a result, for a given and a given , the overall time of our algorithm would increase only linearly with the number of valid cycles and the number of vertices in the graph. In Section 6, we experimentally examine the growth of the running time with respect to , and therefore with , which itself is likely to grow with .
It should be noted that explicit construction of is neither necessary for the correctness of the algorithm, nor does it have a significant effect on its time complexity. By construction, function CYCLE_SEARCH does not visit any vertex whose minimum distance from is or greater. By using instead of , CYCLE_SEARCH avoids testing unnecessary edges at the periphery of the feasible neighborhood of on Line 9. Additionally, explicitly constructing with its own local indices improves the memory access locality of a computer implementation of the algorithm and yields better performance.
6 Experimental results
is the standard deviation of the out-degrees of the vertices.
contains detailed information about each graph. These graphs include social and communication networks, web graphs, product relationships, and professional networks, etc. They range in average density from roughly 1.6 to 18.8. The standard deviation of the out-degrees shows that with the exception of amazon0505, the graphs are fairly nonuniform. This is not surprising, since it has been observed that many social networks are highly skewed graphs and some even exhibit power-law characteristics. For the results reported in this section, a serial implementation of the algorithm in Figures 1–3 was run on a single 3.3 GHz Intel Xeon E5-2667 core. This section contains results from our implementation only because we were unable to find any other practical algorithms or software to solve this problem for comparison on large graphs.
Table 2 shows the number of cycles detected and the running time of our algorithm for = 4, 5, and 6. While our algorithm detects all simple cycles of length less than or equal to , we do not count and report cycles of length one (self-edges) and two (bidirectional edges) in this table. The table shows that the algorithm successfully detected up to hundreds of billions of cycles on large graphs using a single CPU, with only one run (web-BerkStan graph for ) that did not complete in several days and had to be aborted. By default, we sort and renumber the vertices of the graph in increasing order of total degree, as discussed in Section 3.4. For = 5, we also show the running time with this sorting turned off to show the benefit of sorting. The time is worse without sorting for every graph in our test set. In most cases, sorting results in an improvement by a relatively small factor; however, for one graph (wiki-topcats), the difference is so significant that the algorithm may not be practical without sorting.
In Section 5, we derived as the expression for the time complexity of our algorithm for graphs with a uniform out-degree . We tested the graphs in our test suite from the SNAP data set against this expression for growth of running time with and . Table 3 shows the ratio for for these graphs. We used min in the denominator because is meant to represent the maximum number of edges in each in the CYCLE_SEARCH algorithm (Figure 2) for a given and in practice, it cannot exceed the total number of edges in the original graph. Table 3 shows that this ratio decreases as (and therefore, ) increase for all the graphs in our test set. Even though our graphs are quite diverse in origin and density, and the distribution of edges over the vertices is highly nonuniform, they all share this trend. The analytically derived running time bound used a simplifying assumption that the vertices had a uniform degree . Experimentally observing that despite non-uniform degrees, the growth in the running time with and for every test case is slower than that suggested by the analytical bound is therefore an important finding and reinforces the practicality of our algorithm.
7 Concluding remarks and future work
While enumerating simple cycles in directed graphs is a classical computer science problem, in most real applications on very large graphs, only a small fraction of the potentially enormous number of cycles are of interest. There has been relatively little research on discovering interesting cycles and algorithms for these are typically expensive , unless the search space can be substantially narrowed based on the definition of interesting cycles, for example, cycles with obligatory vertices or edges . One class of cycle enumeration and detection problems of practical interest involves finding cycles of a limited length. We devised an efficient and practical algorithm for finding all length-constrained simple cycles in a directed graph.
We analytically and experimentally demonstrated the efficiency of the algorithm and showed that a serial implementation could detect billions of cycles on large graphs with millions on edges in a reasonable amount of time. In the future, it would be interesting to study the results of a parallel implementation on even larger graphs. On the analysis side, it would be interesting to explore the derivation of tighter time bounds without simplifying assumptions about the structure of the graphs.
-  Florian Adriaens, Cigdem Aslay, Tijl De Bie, Aristides Gionis, and Jefrey Lijffijt. Discovering interesting cycles in directed graphs. In CIKM ’19: Proceedings of the 28th ACM International Conference on Information and Knowledge Management, pages 1191–1200, November 2019.
-  N. Alon, R. Yuster, and U. Zwick. Finding and counting given lenth cycles. Algorithmica, 17:209–223, 1997.
-  Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439):509–512, October 1999.
-  Etienne Birmelé et al. Optimal listing of cycles and st-paths in undirected graphs. In Proceedings of the 2013 Annual ACM-SIAM Symposium on Discrete Algorithms, 2013.
-  Anna D. Broido and Aaron Clauset. Scale-free networks are rare. Nature Communications, 10:1017, March 2019.
-  Kevin Buchin, Christian Knauer, Klaus Kriegel, Andre Schulz, and Raimund Seidel. On the number of cycles in planar graphs. In International Computing and Combinatorics Conference, pages 97–107, 07 2007.
-  Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford S. Stein. Introduction to Algorithms (Second Ed.). MIT Press, McGraw-Hill, New York, NY, 2001.
-  Donald B. Johnson. Finding all the elementary circuits of a directed graph. SIAM Journal on Computing, 4(1):77–84, 1975.
-  Steffen Klamt and Axel von Kamp. Computing paths and cycles in biological interaction graphs. BMC Bioinformatics, 10:181, 2009.
-  Jure Leskovec and Andrej Krevl. SNAP Datasets: Stanford large network dataset collection.
-  B. Monien. How to find long paths efficiently. Annals of Discrete Mathematics, 25:239–254, 1985.
-  Robert E. Tarjan. Depth-first search and linear graph algorithms. SIAM Journal on Computing, 1(2):146–160, 1972.
-  Robert E. Tarjan. Enumeration of the elementary circuits of a directed graph. SIAM Journal on Computing, 2(3):211–216, 1973.