{'fraction_non_alphanumeric': 0.04161353888953286, 'fraction_numerical': 0.011069896835516402, 'mean_word_length': 4.163076002393776, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Loop-based Module Detection in Directed Networks

We investigate the problem of finding modules (or clusters, communities) in directed networks. Until now, most articles on this topic have been oriented towards finding full network partitions, where every node belongs to exactly one module. In this paper, we will present a novel random walk based approach for fuzzy partitions of the directed network into modules, such that some nodes do not belong to only one of the modules but to several or to none at all. Despite the fact that the new random walk process is reversible, we will show that it inherits all necessary information about directions and structure of the original network. We accompany our theoretical considerations with numerical examples.

Real-world systems are often modeled as complex networks (or graphs). Understanding the behavior of these systems is thus closely related to investigating the topology and dynamics of the associated networks. One of the most challenging questions in this direction is efficient and accurate identification of network modules, i.e. densely inter-connected subgraphs having sparse connections to the rest of the network. During the last years, different approaches have been proposed to detect modules in networks [1][2][3][4], see [5] for an exhaustive review. However, most of them suffer from the following limitations: (i) they consider only full partitions of a network into modules, i.e. every node is assigned to exactly one module; (ii) they can be applied only to undirected networks. While there is a large variety of methods ranging from pure graph algorithms like mincut via spectral methods to clustering techniques for type (i) of problems, methods for fuzzy partitions of directed networks are still rare. The largest family of such methods seems to be based on Stochastic Blockmodel [6] and Variational Bayesian [7] methods.
The two main results of this paper are: (a) a generalization of the well known modularity function [2] to directed networks which encodes directional information and (b) a new method for finding fuzzy partitions of directed, weighted networks using random walk processes (RW). We consider the network G = (V ; E), where V is the set of nodes and E the set of edges. We further assume that the network is strongly connected, in particular that it has no sinks and sources. Let us denote the weight of edge (xy) ∈ E by K(x, y) and the weighted in-and out-degree of a node by K + (x) = y K(x, y), K − (x) = y K(y, x) respectively. Now, we can define a time-discrete RW process (X n ) n on the network by specifying the following transition matrix P Because of our assumptions on the networks we consider, P is ergodic with unique stationary distribution π. If the network is undirected then K is symmetric and P is reversible.
A widely used approach for finding modules C m ⊂ V in directed networks is to optimize a quantity Q called modularity over all module assignment functions χ m , where χ m (x) = 1 for x ∈ C m and χ m (x) = 0 for x / ∈ C m . The modularity function has been introduced by Newman and Girvan [2,8] for undirected, unweighted networks and various generalizations have been proposed [9][10][11] that extend the definition to directed networks, all of which can be written as Optimization of Q is known to be NP-hard [12], but various fast heuristic algorithms exist [13,14]. We will refer to this family of methods as generalized Newman-Girvan (gNG) approaches. It is obvious that gNG methods suffer from two main drawbacks. First, they consider only full partitions. Second, Q even in its generalized form is actually blind to the directional structure in the network: Q is left unchanged if we replace π x p xy by its symmetrized version π x p s xy = 1 2 (π x p xy + π y p yx ), which is not sensitive to edge directions: π x p s xy can be large if for example only the edge (xy) exists, independent of the length of the shortest path from y back to x. In other words, optimizing Q can produce modules consisting of nodes that communicate only in one direction, which contradicts the intuition that all nodes in a module should be similar.
In this article, we require nodes in a module to communicate in both directions via short paths, and we will make this notion precise by introducing a novel measure of communication I xy between nodes based on cycles. Despite the fact that I xy is symmetric, we will show that it inherits all necessary information about directions and structure of the original network. It will allow us to (a) propose a modified versionQ of the modularity function which is sensitive to directional information and (b) transform G into an undirected network were x and y are connected by an edge with weight I xy if I xy > 0. On the transformed network we can then use any fuzzy clustering method designed for undirected networks.
Cycle decomposition and communication between nodes. We now briefly introduce the theory of cycle de-arXiv:1404.2813v1 [physics.data-an] 10 Apr 2014 compositions for Markov chains as developed in [15] and [16], which we will use for the RW process (1). An ncycle (or n-loop) on G is defined as an ordered sequence [17] of n connected nodes γ = (x 1 , x 2 , . . . , x n ), whose length we denote by |γ| = n. Let C be the collection of simple (i.e. no self-intersections are allowed) cycles on G. Let (X n ) 1≤n≤T be a realization of the Markov chain. We say that (X n ) 1≤n≤T passes through the edge (xy) if ∃n < T such that X n = x and X n+1 = y, and following [15,16] we say that (X n ) 1≤n≤T passes through γ if it passes through all edges of γ in the right order, but not necessarily consecutively. Let N γ T be the number of times γ is passed through up to time T . Then the limit exists almost surely [16] and gives us a uniquely defined probabilistic cycle decomposition, that is a collection Γ = {γ ∈ C|w(γ) > 0} of cycles with positive weights w(γ) such that for every edge (xy) ∈ E the flow decomposition formula holds: where F xy = π x p xy is the flow through (xy) and we write γ ⊃ (xy) if the edge (xy) is in γ. An explicit but impractical formula to calculate the weights w(γ) was given in [16]. Alternatively, and in analogy to how one samples the transition matrix from a realization (X n ) 1≤n≤T , we can also sample them by obtaining the counts N γ T . In [16] such a sampling algorithm was described in detail.
Example: The barbell graph. As an example consider the barbell graph presented in Figure 1, with edge weight K(l 0 , r 0 ) = K(r 0 , l 0 ) = ε < 1, all other edge weights are one. Since every edge belongs to exactly one of the three loops α c = (l 0 , r 0 ), α l = (l 0 , l 1 , . . . , l n−1 ) and α r = (r 0 , r 1 , . . . , r n−1 ), the weights of these loops can be inferred directly from (4): We now define a measure for communication between nodes x, y ∈ V , as experienced by the realization (X n ) 1≤n≤T . To this end, think of (X n ) n as describing a 'postman' delivering parcels. Each time x is visited, a parcel is picked up, which is then delivered when y is reached. Let N T (x → y) be the number of deliveries from x to y in time T -note that x and y need not be adjacent. Then is a measure for the communication between x and y, as experienced by (X n ) n . This quantity is called reaction rate in Transition Path Theory [18]. It is easy to see that k xy = k yx , which does not mean that k xy ignores directional information. In fact, k xy contains information about all possible ways to go from x to y together with information about all possible ways to go from y to x, since the RW must return to x from y before going to y again. Now focus on a single event x → y → x, which is realized by passing through a loop γ. This event gives |γ| counts to y N T (x → y ), one for each y ∈ γ. In other words, the single event x → y → x along γ is counted as one parcel picked up, but |γ| parcels delivered. This 'overcounting' leads to y k xy not being meaningfully normalized.
To arrive at a normalized measure for communication between nodes, we will let the RW pick up the parcel at x and keep it until it returns to x, passing through a loop γ. Then it will be delivered to any node on γ with probability 1 |γ| . Now the numberÑ T (x → y) of deliveries from x to y isÑ Finally, we define the communication intensity I xy as the average rate of parcel deliveries: using (3) and (6), and we arrive at the probabilistic cycle decomposition introduced earlier. Intuitively, I xy is large if there are many cycles connecting x and y, and if they are important (w(γ) large) and short (|γ| small). Thus, our new node communication intensity measure I xy encodes the directional information. Additionally, I xy is normalized as where the last equality follows from (4). This allows us to introduce the loop transition matrixP with components Note thatP is reversible since I xy = I yx and it has the same stationary distribution as P , namely π. Using I xy as a measure for the communication between nodes has provided us with a way to symmetrize P without loosing the directional information.
Generalized modularity function. We now define the modified modularity function which encourages two nodes x,y to belong to the same module if I xy is large, that is if x and y are connected by many short cycles with large weights w(γ). The advantage ofQ compared to Q is that it explicitly requires nodes in the same module to be connected in both directions via short paths. However, optimizingQ will not result in fuzzy partitions. Instead we note that I xy also gives us a transformation of G into an undirected, weighted network G U , where we connect two nodes x, y by an edge with weight I xy if I xy > 0. Now we can use any fuzzy clustering algorithm for undirected networks on G U . Notice that the standard transition matrix (1) on this network coincides withP .
Fuzzy clustering algorithm for directed networks. We now present our new LOop-based Lumping Algorithm (LOLA) for fuzzy clustering of directed weighted networks: (A) Given edge weights K(x, y), construct P according to (1) and generate a realization (X n ) 1≤n≤T of the Markov chain given by P .
(B) Use the sampling algorithm described in [16] to obtain the cycle decomposition Γ with weights w(γ) via (3) and construct G U using (7).
(C) Use any fuzzy clustering method on G U .
In this paper, we will use the Markov State Model (MSM) clustering method [19,20] in (C). MSM clustering identifies module cores C m as the metastable sets of the RW process on G U , which hasP as it's transition matrix. Fuzzy affiliation functions are obtained as q m (x) = P(X t hits C m next|X 0 = x), ∀x ∈ V.
MSM clustering relies on the reversibility ofP [21,22] and can therefore not be applied to directed networks a priori.
Remark 1. If G is very large and/or very metastable, the sampling in step (B) can become too slow. To increase performance, we can replace the sampling algorithm in step (B) by any deterministic algorithm yielding a different cycle decomposition (Γ,w) which will still satisfy (4), but not (3), and will no longer be unique. In [15] a deterministic algorithm was given which can find a decomposition in polynomial time by iteratively reducing the flow F , and the matrixP (Γ) with components is still reversible with stationary distribution π. Remark 2. If the network is undirected, every edge gives rise to a 2-loop. A cycle decomposition using only these loops exists by puttingw((xy)) = π x p xy , and a short calculation showsP (Γ) = 1 2 (I + P ). For MSM clustering, this means that the RW is only slowed down, but will produce the same modules. We now illustrate the method by investigating two small example networks and a real-world network constructed from earthquake data.
The barbell graph cont'd. We can calculate I xy explicitly: x, y ∈ α l w/n x, y ∈ α r εw/2 x = y ∈ α c εw/2 + w/n x = y ∈ α c Clustering this graph with LOLA produces a full partition into two modules C 1 = α l and C 2 = α r , see Figure  1. We will call this partition the reference partition. The scores assigned by Q andQ to it are almost identical: Now we consider what happens if one module is split in two chains of equal size C 1,a and C 1,b . The total change ∆Q in Q resp. ∆Q inQ under this split is One can check that ∆Q > 0 for n ≥ 8 and ∆Q < 0 as long as n > ε. That is, as n grows Q favors a partition into more and more subchains with less then 8 nodes over the reference partition, even though increasing n actually increases the metastability of the reference partition. In contrast, the small chains favored by Q are not metastable at all. On the other handQ always favors the reference partition. Finally, one can check that if directions are ignored, Q behaves in exactly the same way as before, whileQ favors a partitioning into subchains with less then 6 nodes since the nodes in these subchains can now communicate directly. Thus, in this example we can see that Q is not able to detect a qualitative difference between the undirected and directed barbell graph, whileQ is.
False module identification. The next example is a network with 21 nodes, for which LOLA-clustering finds two metastable modules: C 1 colored in blue and C 2 colored in red, see Figure 2a). The rest of the network forms a large transition region consisting of nodes with affiliation less than 0.8 (see Supplementary material for more detailed analysis). If we cluster this network using the gNG algorithm, a third module C 3 appears (green in 2b). However, C 3 is not a metastable module because none of its nodes are connected via short paths in both directions. For example, A and B are connected by a directed edge (AB), but in order to go from B to A, the RW has to pass through the whole network. Consequently, I AB is small, and therefore LOLA-clustering overcomes this problem and is thus improving upon existing methods. Earthquake graph. The last example is a timeseries (X n ) n of seismic events in California from 1952 to 2012, obtained from the SCEC [23]. A weighted, directed network is constructed from this data as follows: Only events with magnitude larger than m c = 2.5 are considered (these are 48669 events). Space is discretized into boxes of ∆l = 0.1 • in the latitude and longitude directions, and the boxes in which earthquakes occurred are the nodes. Then edge weights are defined as K(x, y) = P(X n = x, X n+1 = y), which is estimated from the timeseries by counting successive events. See [24][25][26][27] for a discussion of this and related approaches to construct networks out of seismic data.
The presence of the timeseries data allows us to achieve a significant speed up by constructing the cycle decomposition in step (B) directly from the data. This entirely avoids the construction of P , and no further sampling is needed. For the network in question which has 2175 nodes and 28839 edges, step (B) takes 8 minutes on a standard laptop and reports 7739 cycles. This speedup will always be possible if the network is constructed from timeseries data, which is very common. We give more details in the supplementary material available at [url].
The fuzzy clustering obtained by LOLA is shown in Figure 3, where nodes are colored if they are assigned to a module with affiliation at least 0.8, and otherwise assigned to the transition region and shown grey. In fact 80% of the nodes are assigned to the transition region, but these correspond to only 25% of all events. This illustrates that our fuzzy clustering correctly reflects the uncertainty coming from limited data. A full clustering would have to cluster the grey nodes as well, even though not enough data is available to do so. LOLA-clustering finds 9 modules, all of which correspond to important faults or groups of faults, the largest one containing the San Andreas fault. Although the idea of using networks for analyzing seismic data is not new, to the best of our knowledge, our method is the first that can identify faults based only on the data and without using any additional information.
FIG. 3: Quaternary faults [28] in Southern California and the clustering of the SCEC timeseries found by LOLA. Node size is proportional to the number of events, color indicates the modules found.
Conclusion. In this paper we addressed the problem of module detection in weighted directed networks by constructing a novel measure of communication between nodes. This measure is based on a cycle decomposition of the probability flow which encodes the directional information. Since our measure is symmetric, it allows us to apply clustering methods designed for undirected graphs. We applied our new method on analyzing two toy examples and showed how it overcomes the essential drawbacks of common methods. Finally we demonstrated our novel approach on a real-world directed network coming from a timeseries data.
All RW-based approaches fail when confronted with networks containing sink and source nodes, and our algorithm is no exception, the reason being that the RW is no longer ergodic. A common ad hoc solution is to ergodify the process by including a small teleportation probability τ , this is done e.g. in Google's PageRank algorithm [29]. In future work, we hope to address this