Fuzzy Random Walkers with Second Order Bounds: An Asymmetric Analysis

Abstract

Edge-fuzzy graphs constitute an essential modeling paradigm across a broad spectrum of domains ranging from artificial intelligence to computational neuroscience and social network analysis. Under this model, fundamental graph properties such as edge length and graph diameter become stochastic and as such they are consequently expressed in probabilistic terms. Thus, algorithms for fuzzy graph analysis must rely on non-deterministic design principles. One such principle is Random Walker, which is based on a virtual entity and selects either edges or, like in this case, vertices of a fuzzy graph to visit. This allows the estimation of global graph properties through a long sequence of local decisions, making it a viable strategy candidate for graph processing software relying on native graph databases such as Neo4j. As a concrete example, Chebyshev Walktrap, a heuristic fuzzy community discovery algorithm relying on second order statistics and on the teleportation of the Random Walker, is proposed and its performance, expressed in terms of community coherence and number of vertex visits, is compared to the previously proposed algorithms of Markov Walktrap, Fuzzy Walktrap, and Fuzzy Newman–Girvan. In order to facilitate this comparison, a metric based on the asymmetric metrics of Tversky index and Kullback–Leibler divergence is used.

1. Introduction

The Random Walker principle is the algorithmic cornerstone for building a number of heuristics for large graphs, namely for those with the fundamental property that neither their vertex nor their edge set fits in main memory. Such heuristics are efficient in terms either of computation time or memory requirements or often both. Under this principle, a virtual entity usually named the Random Walker visits the vertices. Within the scope of this article, the probabilistic strategy followed by the Random Walker to decide which vertex will visit next is of paramount importance, although, depending on the problem under study, other properties of the Random Walker may be of interest.

Virtual or ideal entities play an important role in science and engineering, mainly as a means to prove a theorem, to establish ideal performance limits, and to provide grounds for rejecting a conjecture based on a reductio ad absurdum methodology. Consider, for instance, the particle sorting demon of Maxwell [1,2] with its connections to algorithmic information theory and the steam engine of Heron of Alexandria [3]. In addition, the Random Walker principle itself has been applied to a number of graph analytics such as vertex similarity [4] and graph cuts [5] as well as to image processing [6].

A central finding of the theory of scale-free graphs, corroborated by significant evidence from a broad range of fields such as genomics, computational neuroscience, combinatorics, and social network analysis, states that they strongly tend to exhibit modularity. Namely, scale free graphs are recursively composed of vertex communities. The latter is a crucial factor for simultaneously achieving scalability, low-delay local information propagation, and structural robustness [7,8]. For instance, a large online social media group for Roman history afficionados can be composed of overlapping but to a great extent distinct specialized communities about Roman politics, law, and military, whereas the latter may be further subdivided into smaller communities regarding weaponry, tactics, legion standards, key battles, distinguished commanders, and the Praetorian guard.

Even though knowledge about communities offers a deep insight to graph structure, locating them is intractable. Consequently, numerous community discovery algorithms have been developed. Two of the most prominent ones are the Newman–Girvan [8] and the Walktrap [9] algorithms. The former is deterministic and is based on local edge density, while the latter is heuristic and relies on an edge crossing Random Walker. In [10], variants of both algorithms for edge-fuzzy graphs were proposed. Furthermore, in [11], a random walk approach for community detection in complex networks is introduced.

This article extends the work of [12]. Initially, the implementation of the Markov Walktrap algorithm, which is an extended version of Fuzzy Walktrap, was firstly proposed in [10]. Fuzzy Walktrap in turn has been based on deterministic Walktrap algorithm [9]. Furthermore, a method for algorithmically assessing the performance of Markov Walktrap relative to Fuzzy Walktrap and Fuzzy Newman–Girvan, being another fuzzy community detection technique [10], are thoroughly presented. Finally, one very important aspect to be mentioned is the expression of Markov Walktrap in Cypher, which composes the main query language for Neo4j [13].

The primary contribution of this article is the implementation over Neo4j of Chebyshev Walktrap, a community discovery algorithm designed for edge-fuzzy graphs, a class of fuzzy graphs [14] used among others in [10], which is based on the Random Walker algorithmic principle. Additionally, Chebyshev Walktrap relies on the competitive factors of second order statistics though the Chebyshev inequality and on an optional relocation capability in order to bound unnecessarily costly walks and, thus, remaining inside a community and being trapped for too long within the boundaries of a community, respectively. The relocation aspect was also backported to the Markov Walktrap algorithm first proposed in [15]. The effect of relocation on the community coherence was evaluated based on the asymmetric Tversky index using the Fuzzy Newman–Girvan algorithm from [15] as baseline, while its effect on the output distribution was assessed with the asymmetric Kullback–Leibler divergence.

The remainder of this article is structured as follows. Scientific literature regarding community discovery is reviewed in Section 2. The fuzzy graph model is outlined in Section 3 while Markov Walktrap and Fuzzy Walktrap algorithms are presented in Section 4. Experimental results are discussed in Section 5. Finally, the main findings of this paper as well as future research directions are discussed in Section 6.

Table 1. Paper notation.

Symbol	Meaning
$\stackrel{▵}{=}$	Definition or equality by definition
⊗	Kronecker tensor product
$\left\{s_1,\dots ,s_n\right\}$	Set with elements $s_1,s_2,\dots ,s_n$
$\left|·\right|$	Set cardinality or path length (depending on the context)
$\left(e_1,\dots ,e_m\right)$	Path comprised of edges $e_1,\dots ,e_m$
$K_n$	Complete graph with n vertices and $\left(\genfrac{}{}{0pt}{}{n}{2}\right)$ edges
$\mathrm{E}\left[X\right]$	Mean value of random variable X
$Var\left[X\right]$	Variance of random variable X
${\tau }_{T,V}$	Tanimoto similarity coefficient for sets T and V
${\nu }_{T,V}$	Asymmetric Tversky index for sets T and V
$S_1\setminus S_2$	Asymmetric set difference $S_1$ minus $S_2$
$\tilde{S}$	Fuzzy set S
$〈s_k〉$	Sequence of elements $s_k$
$\mathcal{H}\left(s_1,\dots ,s_n\right)$	Harmonic mean of elements $s_1,\dots ,s_n$
$\mathcal{H}\left(s_1,\dots ,s_n;{\tau }_0\right)$	Thresholded or effective harmonic mean of $s_1,\dots ,s_n$
${\underline{\mathbf{1}}}_n$	$n\times 1$ vector with ones
${\underline{\mathbf{e}}}_n^k$	$n\times 1$ zero vector with a single one at the k-th entry
$f^{\left(n\right)}\left(x\right)$	n-th order derivative of $f\left(x\right)$
$〈p\left|\right|q〉$	Kullback–Leibler divergence between distributions p and q

summarizes the symbols used in this paper.

2. Related Work

From a system perspective, graph analytics play a crucial role in data mining systems such as Google Pregel or in massive machine learning frameworks such as GraphLab (https://turi.com/). Moreover, recently, there is strong interest for scalable, production grade graph databases [13,16] such as BrightStar (https://brightstardb.com), Neo4j (www.neo4j.com), Titan (https://github.com/thinkaurelius/titan), Sparksee (www.sparsity-technologies.com) and GraphDB (www.ontotext.com).

Traditionally, from an algorithmic viewpoint, analytics include structural [17,18] and spectral [19,20] partitioning, where a graph is split according to some functional constraints such as flow or edge density. Efficient information diffusion in large graphs is also of interest [21,22], especially for online political campaigns and digital marketing. The Random Walker principle has been also applied to two other important metrics, namely vertex similarity [4] and heuristic minimum cuts [5]. Both metrics can also be treated deterministically, especially in the context of social network analysis [23,24]. Community structure discovery provides insight to the inner workings of a particular graph [7,8,17], while metrics such as those in [25] control the discovery process quality. Persistent graphs can be instrumental in designing rollback capabilities in graph databases [26].

Among the numerous applications of graphs or linked data, one can find Web searching and ranking [27] with established algorithms such as PageRank [28] and HITS [29]. Bibliometric and scientometric data analysis [30] can boost collaboration between researchers, while image segmentation [6] is central to computer vision and robotics. Social network analysis has greatly benefited from structural [31] or functional [32,33] community detection algorithms. Additionally, influence and perceived social status in online social media have been tied to the participation in communities [34]. Message diffusion within a social graph is studied [34,35], while [36,37,38] deal with emotional modeling with respect to user influence [36]. Finally, random walkers have served as models for the propagation of computer viruses both in single systems and in networks, including LANs and the Internet [39,40]. Within this context, the strategy or the mix of strategies followed by the random walker is of paramount importance as it affects the entity and number of resources susceptible to infection.

First and second order statistics are used across a number of fields. In [41], a channel estimation methodology based on first order statistics is proposed. Methods for blind source separation using second order statistics include [42,43]. A comprehensive approach about the applications of higher order methods is given in [44] signal processing and in [45] for biomedical engineering. In [46], a third order method was presented for adaptively scheduling biosignal processing applications at the operating system level. Independent Component Analysis (ICA), a powerful signal processing technique, is based on higher order spectra [47]. Among the multitude of ICA applications is source separation in EEG waveforms [48].

3. Edge-Fuzzy Graphs

3.1. Definitions

Within the scope of this paper, the edge-fuzzy graphs are probabilistic and combinatorial hybrids comprised of a fixed set of vertices V and a fuzzy set of edges $\tilde{E}$. Formally,

Definition 1.

A homogeneous edge-fuzzy graph is the ordered triplet

$$G\stackrel{▵}{=}\left(V,\tilde{E},h\right),$$

(1)

where $V=\left\{v_k\right\}$ is the set of vertices, $\tilde{E}=\left\{e_k\right\}\subseteq V\times V$ is the fuzzy set of edges, and h is the edge membership function $h:E\to \left(0,1\right]$, which quantifies the degree of participation of each $e_k$ to G [10]. Moreover:

• the vertex set V is fixed, namely they belong to G with probability one,
• the distribution h is the same for each edge,
• the existence probability of $e_k$ is drawn independently for each edge.

A subtle point is that h does not affect the structural properties of the graph in the sense that no edges are added or deleted, except when for a particular $e_k$ holds that $h\left(e_k\right)=0$. If h is continuous, the probability of this ocurring is zero. However, if h is discrete, then depending on h a potentially non-negligible portion of the edges may be deleted. In this article, h was chosen so that only at most an exponentially small proportion of the edges would be assigned to a zero weight. Consequently, the underlying graph preserved its original structure along with any associated connectivity patterns. Otherwise, if a considerable fraction of edges were to be deleted, then the resulting graph would behave more like an Erdös-Rényi graph. The latter are known to be easily constructed by randomly sampling a graph space or, equivalently, the edges of $K_n$ but their properties deviate in a significant way from those of real world, large graphs.

Observe that there is no fuzziness whatsoever regarding vertices as by definition they always exist with probability one. In scientific literature, the existence of fuzzy graph classes is prominent. Vertices are fuzzy as well as their fuzziness interacts with that of the edges, mostly by long product chains. Such graphs are beyond the scope of this article.

Definition 2.

Under the fuzzy graph model of Defintion 1, the cost $\delta \left(e_k\right)$ of traversing $e_k$ is

$$\delta \left(e_k\right)\stackrel{▵}{=}\frac{1}{h\left(e_k\right)}\in \left[1,+\infty \right),h\left(e_k\right)\in \left(0,1\right],$$

(2)

which expresses the intuitive requirement that edges which are less likely to belong to the graph are also harder to cross.

Depending on the application, $\delta$ may well be connected through another non-linear transform to h as long as edges with high $h\left(e_k\right)$ are easy to cross and edges with low h are difficult to cross such as

$$\begin{array}{c}{\delta }^{\prime }\left(e_k\right)\stackrel{▵}{=}\frac{1}{1+h^2\left(e_k\right)}\in \left[\frac{1}{2},1\right],\hfill \\ {\delta }^{″}\left(e_k\right)\stackrel{▵}{=}ln\left(\frac{1}{h\left(e_k\right)}\right)=-lnh\left(e_k\right)\in \left[0,+\infty \right).\hfill \end{array}$$

(3)

Definition 3.

The cost $\Delta \left(p_j\right)$ of a fuzzy path $p_j=\left(e_1,\dots ,e_m\right)$ is the sum of the cost of its individual edges

$$\Delta \left(p_j\right)\stackrel{▵}{=}\sum _{e_k\in p_j}\delta \left(e_k\right)=\sum _{k=1}^m\frac{1}{h\left(e_k\right)}=\frac{m}{\mathcal{H}\left(h\left(e_1\right),\dots ,h\left(e_m\right)\right)},h\left(e_k\right)\ne 0,1\le k\le m,$$

(4)

where $\mathcal{H}\left(h\left(e_1\right),\dots ,h\left(e_m\right)\right)$ is the harmonic mean of $h\left(e_k\right)$ defined as

$$\mathcal{H}\left(s_1\dots s_n\right)\stackrel{▵}{=}\frac{1}{\frac{1}{n}{\sum }_{k=1}^n\left(\frac{1}{s_k}\right)}=\frac{n}{\frac{1}{s_1}+\dots +\frac{1}{s_n}},s_k\ne 0,1\le k\le n.$$

(5)

By construction $\Delta \left(p_j\right)$ is bounded as follows

$$\frac{1}{{max}_{1\le k\le m}\left\{h\left(e_k\right)\right\}}\le \frac{\Delta \left(p_j\right)}{\left|p_j\right|}\le \frac{1}{{min}_{1\le k\le m}\left\{h\left(e_k\right)\right\}}.$$

(6)

Whether the above bounds are loose depends on the variability of the actual values $h\left(e_k\right)$. That is, if the latter are drawn from a distribution which favors extreme values, $\Delta \left(p_j\right)$ will tend to be close to these bounds. It should also be noted that the variance $\Delta \left(p_j\right)$ is strongly dependent on that of $h\left(e_k\right)$. Moreover, $\Delta \left(p_j\right)$ is prone to outliers, which might lead to an unrealistically high average fuzzy path length. This can be remedied by taking into account the variance of the fuzzy path length. Finally, a low $\Delta \left(p_j\right)$ tends to contain almost exclusively low edge costs $\delta \left(e_k\right)$ or, equivalently, edges with high probability of belonging to the graph, an argument which agrees with the weak law of large numbers. In turn, this suggests the intuitive corollary that low cost paths are comprised almost exclusively of edges that are less likely to be fuzzy. This corollary can be used in order to design efficient hybrid probabilistic and combinatorial algorithms based on dynamic programming for finding and enumerating low cost paths akin to the way a similar observation has led to the development of shortest paths relying on dynamic programming in deterministic graphs.

From a probabilistic viewpoint, the sum ${\sum }_{k=1}^m\delta \left(e_k\right)$ is interesting by itself as a finite but possibly large sum of inverse random variables. Notice that the central limit theorem may not be applied in such a setting since the variance of the $h\left(e_k\right)$ might be infinite. If this is not the case, different bounds can be computed depending on the distribution of $h\left(e_k\right)$ such as the abovementioned central limit theorem, a Poisson bound, a power law bound, or finally approximations based on Markov or Chebyshev inequalities, the Chernoff bound, or on the Gnedenko extreme value theorem. Notice that the effect of a single edge which is exceedingly difficult to cross can be instrumental in shaping graph communities.

The numerical properties of the above sum are also of interest. As values of possibly uneven orders of magnitude may be added, catastrophic cancellation may occur resulting in the loss of meaningful information. This might happen if the summation is executed in an order left to the implementation. On the other hand, the summation order dictated by the Priest algorithm [49] results in the least possible loss of significant decimal digits by adding only numbers of comparable magnitude. Another option would be to substitute the harmonic mean with its thresholded counterpart

$$\mathcal{H}\left(s_1,\dots ,s_m;{\tau }_0\right)\stackrel{▵}{=}\frac{m}{{\sum }_1^m\frac{1}{max\left\{h\left(e_k\right),{\tau }_0\right\}}}=\frac{m}{\frac{1}{max\left\{h\left(e_1\right),{\tau }_0\right\}}+\dots +\frac{1}{max\left\{h\left(e_m\right),{\tau }_0\right\}}}.$$

(7)

For other uses of the thresholded harmonic and geometric means, see [50], while, for the effect of finite precision arithmetic to long biosignals, see [51].

An alternative for long paths would be to substitute, under certain conditions, the finite sum with an appropriate integral. Assuming with no loss of generality that $h\left(e_1\right)={min}_{1\le k\le m}\left\{h\left(e_k\right)\right\}\ne 0$ and $h\left(e_m\right)={max}_{1\le k\le m}\left\{h\left(e_k\right)\right\}\ne 0$, then

$$\Delta \left(p_j\right)\stackrel{▵}{=}\sum _{k=1}^m\frac{1}{h\left(e_k\right)}\approx {\rho }_0+{\int }_{h\left(e_1\right)}^{h\left(e_m\right)}\frac{1}{x}dx=lnh\left(e_m\right)-lnh\left(e_1\right)+{\rho }_0=ln\left(\frac{h\left(e_m\right)}{h\left(e_1\right)}\right)+{\rho }_0,$$

(8)

where ${\rho }_0$ is an optional correction factor. A finer approach requiring more probabilistic information about the longer paths of a given graph would be to partition such a path $p_j$ so that

$$\Delta \left(p_j\right)\stackrel{▵}{=}\sum _{k=1}^m\frac{1}{h\left(e_k\right)}\approx \sum _{i=1}^n\left(\frac{lnh\left(e_{u_i}\right)}{lnh\left(e_{{\ell }_i}\right)}+{\rho }_i\right)=\sum _{i=1}^n\left(\frac{lnh\left(e_{u_i}\right)}{lnh\left(e_{{\ell }_i}\right)}\right)+\sum _{i=1}^n{\rho }_i=\sum _{i=1}^n\left(\frac{lnh\left(e_{u_i}\right)}{lnh\left(e_{{\ell }_i}\right)}\right)+{\rho }_0^{\prime }.$$

(9)

Selecting n and forming the sets ${\left\{h\left(e_{u_i}\right)\right\}}_{i=1}^n$, ${\left\{h\left(e_{{\ell }_i}\right)\right\}}_{i=1}^n$, and ${\left\{{\rho }_i\right\}}_{i=1}^n$ is not a trivial task. Instead, choosing such an approach might be a viable solution only for certain combinations of h and $\left|p_j\right|$. Techniques for estimating the variability as well as the cardinality of large sets such as [52] can be useful while pursuing this approach.

It should be emphasized that the class of fuzzy graphs of Definition 1 can be well considered as a typical example of higher order data. This is attributed to the inherently distributed way information is stored in a graph, in this particular case as edge existence probabilities. In order for meaningful information regarding path costs to be mined, a non-negligible fraction of edges must be crossed and, thus, the interplay of a number of edges must be considered.

3.2. Reciprocal Random Variables

Because of the Definitions 2 and 3 for $\delta \left(e_k\right)$ and $\Delta \left(p_j\right)$, respectively, the properties of an inverse random variable gain more interest. The following definition is straightforward.

Definition 4.

The inverse distribution of a mass distribution function of a random variable X is defined as the mass distribution function of $\frac{1}{X}$ [53].

Property 1.

In the continuous case, the distributions of X and $\frac{1}{X}$ are linked as

$$f_{\frac{1}{X}}\left(y\right)=\frac{1}{y^2}{f}_X\left(\frac{1}{y}\right),y\ne 0.$$

(10)

Proof. 

The cumulative distribution of $\frac{1}{X}$ is defined as

$$F_{\frac{1}{X}}\left(y\right)\stackrel{▵}{=}prob\left\{\frac{1}{X}\le y\right\}=prob\left\{X\ge \frac{1}{y}\right\}=1-prob\left\{X<\frac{1}{y}\right\}=1-F_X\left(\frac{1}{y}\right).$$

(11)

By differentiating the last relationship, the stated result follows. ☐

For instance, if X is the continuous uniform random variable in $\left[{\alpha }_1,{\alpha }_2\right]$, where ${\alpha }_1,{\alpha }_2\ne 0$, then the distribution of $\frac{1}{X}$ is

$$f_{\frac{1}{X}}\left(y\right)=\frac{1}{y^2\left({\alpha }_2-{\alpha }_1\right)},y\in \left[\frac{1}{{\alpha }_2},\frac{1}{{\alpha }_1}\right].$$

(12)

Despite its simplicity, in certain scenaria, relationship (10) cannot be used. For instance, only the first moments of X may be known or $y=0$ might be a legitimate value, in which case there is a singularity in the inversion of $f_X\left(x\right)$. Instead, bounds are sought for the first moments of $\frac{1}{X}$, which makes more sense from a programming viewpoint in the case of large graphs.

Jensen inequality provides a straightforward way to bound the expected value $\mathrm{E}\left[\frac{1}{X}\right]$ by using the expected value $\mathrm{E}\left[X\right]$ of the non-zero random variable X.

Theorem 1.

(Jensen inequality) For any random variable and any convex function $g\left(x\right)$ provided that both the domains of X and $\mathrm{E}\left[X\right]$ are subsets of the domain of $g\left(·\right)$

$$g\left(\mathrm{E}\left[X\right]\right)\le \mathrm{E}\left[g\left(X\right)\right].$$

(13)

Corollary 1.

The mean value of the strictly positive random variable $\frac{1}{X}$ has a lower bound of

$$\frac{1}{\mathrm{E}\left[X\right]}\le \mathrm{E}\left[\frac{1}{X}\right],X\ne 0.$$

(14)

Property 2.

Function $g\left(x\right)=\frac{1}{x}$ is convex when $x>0$.

Proof. 

Notice that the second derivative $g^{\left(2\right)}\left(x\right)=\frac{2}{x^3}$ is positive when x is positive. An alternative way to prove this claim is to apply the standard convexity definition. For every ${\alpha }_0\in \left[0,1\right]$, $x_1>0$, and $x_2>0$

$$\begin{array}{c}\hfill g\left({\alpha }_0{x}_1+\left(1-{\alpha }_0\right)x_2\right)\le {\alpha }_{0}g\left(x_1\right)+\left(1-{\alpha }_0\right)g\left(x_2\right)\Rightarrow \\ \hfill \frac{1}{{\alpha }_0{x}_1+\left(1-{\alpha }_0\right)x_2}\le \frac{{\alpha }_0}{x_1}+\frac{1-{\alpha }_0}{x_2}\iff \\ \hfill {\alpha }_0{x}_2\left({\alpha }_0{x}_1+\left(1-{\alpha }_0\right)x_2\right)+\left(1-{\alpha }_0\right)x_1\left({\alpha }_0{x}_1+\left(1-{\alpha }_0\right)x_2\right)-x_1{x}_2\ge 0\iff \\ \hfill {\alpha }_0\left(1-{\alpha }_0\right){\left(x_1-x_2\right)}^2\ge 0.\end{array}$$

(15)

☐

In order to derive realistic upper bounds for the path lengths of a given fuzzy graph, certain probabilistic inequalities can be employed. The first is Markov inequality, which establishes a first order bound for the probability of X taking very large values by stating that

Theorem 2.

(Markov inequality) The probability of a strictly positive random variable X exceeding ${\gamma }_0$ is bounded by

$$prob\left\{X\ge {\gamma }_0\right\}\le \frac{\mathrm{E}\left[X\right]}{{\gamma }_0},{\gamma }_0>0,X>0.$$

(16)

Second order bounds can be derived by the Chebyshev inequality. The latter provides tighter bounds while lifting the positivity assumption.

Theorem 3.

(Chebyshev inequality) The probability of an arbitrary random variable X exceeding its expected value by a certain fraction ${\gamma }_0$ of its standard deviation as

$$prob\left\{\left|X-\mathrm{E}\left[X\right]\right|\ge {\gamma }_0\sqrt{Var\left[X\right]}\right\}\le \frac{1}{{\gamma }_0^2},{\gamma }_0>0.$$

(17)

The Chebyshev inequality is generic enough to be applied in a number of scenaria, including those in the present article. Still, it should be noted that other techniques may provide sharper bounds in certain cases. For instance, when X is normally distributed, then specialized methods exist for evaluating the integral under its tail.

Estimating the variance of a transformed random variable can be done through the delta method.

Theorem 4.

(Delta method) Let X be a random variable whose expected value $\mathrm{E}\left[X\right]$ and variance $Var\left[X\right]$ are known. For an analytic $g\left(x\right)$, the variance of $g\left(X\right)$ can be estimated as

$$Var\left[g\left(X\right)\right]\approx {\left(g^{\left(1\right)}\left(\mathrm{E}\left[X\right]\right)\right)}^{2}Var\left[X\right].$$

(18)

Proof. 

The first order approximation of the Taylor expansion of $g\left(·\right)$ around $\mathrm{E}\left[X\right]$ is

$$g\left(X\right)=\sum _{k=0}^{+\infty }{g}^{\left(k\right)}\left(\mathrm{E}\left[X\right]\right)\frac{X^k}{k!}\approx g\left(\mathrm{E}\left[X\right]\right)+g^{\left(1\right)}\left(\mathrm{E}\left[X\right]\right)\left(X-\mathrm{E}\left[X\right]\right).$$

(19)

Taking the variance of both sides along with the identities,

$$\begin{array}{cc}\hfill Var\left[{\alpha }_0+{\alpha }_{1}X\right]& ={\alpha }_1^{2}Var\left[X\right],\hfill \\ \hfill Var\left[X-\mathrm{E}\left[X\right]\right]& =Var\left[X\right],\hfill \end{array}$$

(20)

yields the stated result. ☐

Corollary 2.

For $g\left(x\right)=\frac{1}{x}$, the delta method yields

$$Var\left[\frac{1}{X}\right]\approx \frac{Var\left[X\right]}{\mathrm{E}{\left[X\right]}^4}=\frac{\mathrm{E}\left[X^2\right]-\mathrm{E}{\left[X\right]}^2}{\mathrm{E}{\left[X\right]}^4}=\frac{\mathrm{E}\left[X^2\right]}{\mathrm{E}{\left[X\right]}^4}-\frac{1}{\mathrm{E}{\left[X\right]}^2},\mathrm{E}\left[X\right]\ne 0.$$

(21)

The Markov and Chebyshev are but two of the probabilistic inequalities collectively known as concentration inequalities, the latter bound the deviation of a random variable or a sequence of random variables from a known value. Other such inequalities include the Talagrand, the Efron–Stein, and the Dvoretzky–Kiefer–Wolfowitz inequalities.

4. Family of Walktrap Heuristics

4.1. Deterministic Walktrap

The original Walktrap algorithm [9] simulates an edge crossing random walker in order to estimate the stationary distribution of a homogeneous Markov chain. The walker can commence from any vertex and cross edges by randomly selecting destination vertices, systematically moving to vertices with high edge density as vertices are selected with probability proportional to their degree. Since vertices can be visited an arbitrary number of times, unlike algorithms like BFS and DFS, eventually some patterns in the vertex visiting sequence will emerge. As a community is from a structural perspective, essentially a locally dense graph segment, the walker is more likely to move along vertices belonging to the same community for a large time interval before moving to another community. Thus, analysis of the vertex sequence generated by the random walker can reveal the underlying graph community structure. The Walktrap algorithm is outlined in Algorithm 1.

Algorithm 1: Deterministic Walktrap

Require: graph $G(V,E)$, termination criterion ${\tau }_0$ Ensure: vertex pair sequence $〈s_k,s_k^{*}〉$ is generated 1: pick a random vertex $v$ 2: repeat 3:  pick a neighboring vertex $v^{*}$ with probability proportional to its degree as in (23) 4:  store current vertex in v and move the walker to the new vertex $v^{*}$  and $〈s_{k+1},s_{k+1}^{*}〉$ $\leftarrow$ $(v,v^{*})$ 5: until ${\tau }_0$ is true 6: return vertex sequence $〈s_k,s_k^{*}〉$

The degree of any neighboring vertex can be determined by the graph adjacency matrix $\mathbf{A}$ defined as

$$\mathbf{A}\left[i,j\right]\stackrel{▵}{=}\left\{\begin{array}{cc}1,\hfill & i=j\vee \left(i,j\right)\in E\hfill \\ 0,\hfill & i\ne j\wedge \left(i,j\right)\notin E\hfill \end{array}\right.\in {\left\{0,1\right\}}^{\left|V\right|\times \left|V\right|}.$$

(22)

Specifically, the degree of $v_k$ is the sum of the k-th column of $\mathbf{A}$. The probability that from vertex $v_p$ a neighbor $v_q$ is selected at the next step is directly proportional to

$$\frac{deg\left(v_q\right)}{{\sum }_{\left(v_p,v_j\right)\in E}deg\left(v_j\right)}=\frac{{\underline{\mathbf{1}}}_{\left|V\right|}^T\mathbf{A}{\underline{\mathbf{e}}}_{\left|V\right|}^q}{{\sum }_{\left(v_p,v_j\right)\in E}{\underline{\mathbf{1}}}_{\left|V\right|}^T\mathbf{A}{\underline{\mathbf{e}}}_{\left|V\right|}^j}.$$

(23)

In contrast to many graph algorithms, each vertex may be visited more than once. In fact, vertices must be visited many times in order for meaningful patterns to emerge regarding community structure. Typically, for deterministic graphs, a constant number of visits per vertex may suffice resulting in a total of $\mathrm{O}\left(\left|V\right|\right)$ visits, though techniques exploiting the self-similarity nature of large, scale free graphs may yield somewhat lower bounds of $\mathrm{O}\left({log}^{1+ϵ}\left|V\right|\right)$, $ϵ>0$. For fuzzy graphs, the linear bound is as of yet unknown as to whether it can be improved.

In a distributed setting such as Hadoop, the Deterministic Walktrap can be scaled up since graph segments can be distributed to the nodes. The map part will be the parallel random walkers crossing edges. If such a walker must cross a segment, it can either bounce back or be transferred to the appropriate node. The reduce part will be the frequency count of a large number of vertex pairs.

Once the random walker has finished crossing the graph, the communities are discovered by means of hierarchical clustering using the frequency of pairs $\left(s,s^{*}\right)$ as weights. It should be noted that other methods such as Hidden Markov Models and text mining techniques dealing with missing values [54] may be used for discerning community patterns in the sequence $〈s_k,s_k^{*}〉$.

4.2. Fuzzy Walktrap

The Fuzzy Walktrap algorithm has been proposed and analyzed in [10]. Similarly to Algorithm 2, given a vertex $v_p$, each neighbor $v_q$ is a candidate for being visited by the Random Walker with probability proportional to its probability of belonging to the fuzzy graph, namely proportional to

$$\frac{h\left(\left(v_p,v_q\right)\right)}{{\sum }_{\left(v_p,v_j\right)\in \tilde{E}}h\left(\left(v_p,v_j\right)\right)}=\frac{{\underline{\mathbf{1}}}_{\left|V\right|}^T{\mathbf{A}}^F{\underline{\mathbf{e}}}_{\left|V\right|}^q}{{\sum }_{\left(v_p,v_j\right)\in E}{\underline{\mathbf{1}}}_{\left|V\right|}^T{\mathbf{A}}^F{\underline{\mathbf{e}}}_{\left|V\right|}^j},$$

(24)

where the fuzzy adjacency matrix ${\mathbf{A}}^F$ defined as

$${\mathbf{A}}^F\left[i,j\right]\stackrel{▵}{=}\left\{\begin{array}{cc}1,\hfill & i=j,\hfill \\ h\left(e_{ij}\right),\hfill & e_{ij}=\left(v_i,v_j\right)\in \tilde{E}\hfill \\ 0,\hfill & \left(v_i,v_j\right)\notin \tilde{E}.\hfill \end{array}\right.\in {\left[0,1\right]}^{\left|V\right|\times \left|V\right|},$$

(25)

Fuzzy Walktrap is outlined in Algorithm 2.

Algorithm 2: Fuzzy Walktrap

Require: fuzzy graph $G(V,\tilde{E},h)$, termination criterion ${\tau }_0$ Ensure: vertex pair sequence $〈s_k,s_k^{*}〉$ is generated 1: pick a random vertex $v$ 2: repeat 3:  pick a neighboring vertex $v^{*}$ with probability proportional to its cost as in (24) 4:  store current vertex in v and move the walker to the new vertex $v^{*}$  and $〈s_{k+1},s_{k+1}^{*}〉$ $\leftarrow$ $(v,v^{*})$ 5: until ${\tau }_0$ is true 6: return $〈s_k,s_k^{*}〉$

4.3. Markov Walktrap and Chebyshev Walktrap

Both the Markov Walktrap, proposed in [10], and Chebyshev Walktrap, introduced in this article, algorithms improve Fuzzy Walktrap in two ways. The first is that during the walking phase the Random Walker has two optional safeguards against being confined inside a community for too long. Both of these safeguards are common for Markov Walktrap and Chebyshev Walktrap. The second improvement is that, during the clustering phase, two communities may not merge if the path lengths within the resulting community exceed a certain threshold. The latter is based on first order statistics for the Markov Walktrap and on second order statistics for the Chebyshev Walktrap.

Regarding the control of community merge, for community $V_k$, the mean path cost ${\pi }_k$ is the sum of the individual edge costs

$${\pi }_k\stackrel{▵}{=}\frac{1}{\left|V_k\right|}\sum _{e_j\in V_k}\delta \left(e_j\right),$$

(26)

and, therefore, is also a random variable. By linearity of expectation,

$$\mathrm{E}\left[{\pi }_k\right]=\mathrm{E}\left[\frac{1}{\left|V_k\right|}\sum _{e_j\in V_k}\delta \left(e_j\right)\right]=\frac{\left|E_k\right|}{\left|V_k\right|}\mathrm{E}\left[\delta \left(e_j\right)\right].$$

(27)

Moreover, as $\delta \left(e_j\right)$ are independent,

$$Var\left[{\pi }_k\right]=Var\left[\frac{1}{\left|V_k\right|}\sum _{e_j\in V_k}\delta \left(e_j\right)\right]=\frac{\left|E_k\right|}{{\left|V_k\right|}^2}Var\left[\delta \left(e_j\right)\right].$$

(28)

In the general case, the distribution of $\delta \left(e_j\right)$ is unknown, for specific choices of h, it can be computed or estimated. Alternatively, since ${\pi }_k$ is a sum of random variables, its distribution may be known for certain special cases. For instance, the sum of independent Poisson random variable is another Poisson random variable. In addition, the sum of independent binomial random variables is also a binomial random variable. Finally, the sum of a large number of independent random variables with finite variance is a normal random variable according to the Central Limit Theorem. Nonetheless, in this article, no such assumptions were made and the expected value and the variance of ${\pi }_k$ were approximated by the Jensen inequality and the delta method, respectively, in Equations (13) and (18).

Since ${\pi }_k$ is by construction positive, the Markov inequality can be applied. Therefore,

$$prob\left\{{\pi }_k\ge {\gamma }_0\right\}\le \frac{\mathrm{E}\left[{\pi }_k\right]}{{\gamma }_0}$$

(29)

or, equivalently,

$$prob\left\{{\pi }_k\ge \frac{\mathrm{E}\left[{\pi }_k\right]}{{\gamma }_0}\right\}\le {\gamma }_0.$$

(30)

If, for a threshold ${\alpha }_0\in \left(0,1\right)$${\pi }_k$ exceeds ${\alpha }_0{\gamma }_0$, then that community is excluded from merging for an iteration provided it has at least ${\xi }_0$ vertices. This is a first order probabilistic safeguard preventing almost formed communities from losing their coherence.

On similar grounds, a second order such safeguard can be built on the Chebyshev inequality

$$prob\left\{\left|{\pi }_k-\mathrm{E}\left[{\pi }_k\right]\right|\ge {\gamma }_1\sqrt{Var\left[{\pi }_k\right]}\right\}\le \frac{1}{{\gamma }_1^2}.$$

(31)

4.4. Escape Strategies

Although the purpose of the Random Walker is to discover communities by repeatedly visiting neighboring vertices and crossing low cost edges, it is possible to be trapped inside a community if the latter is connected only through very high cost edges from the remaining graph. To this end, the Random Walker has the option as in [15] to reverse its strategy and select neighboring vertices with probability inversely proportional to their probability of existence if a random flag is triggered. The latter was implemented as a Bernoulli random variable with success probability $q_0$, which can be set to zero if so desired in order to disable the weight inversion strategy. Recommended values are typically $\mathrm{O}\left({\left|V\right|}^{-ϵ}\right)$, $ϵ\ge 2$. Therefore, as in [15], when weight inversion is enabled, the distance d between communities implicitly depends on terms of the form

$$d\propto \frac{{\prod }_{e_k\in \tilde{E}}\delta \left(e_k\right)}{{\prod }_{e_j\in \tilde{E}}\delta \left(e_j\right)},$$

(32)

instead of terms of the form

$$d\propto \prod _{e_k\in \tilde{E}}\delta \left(e_k\right).$$

(33)

Alternatively, a probabilistically triggered restart of the Random Walker akin to the PageRank teleportation [28], the random mutation operator in a genetic algorithm [55,56], or the restart strategy in GMRES iterative solver for linear systems [57,58] was also considered. The relocation probability $q_1$ has a Bernoulli distribution and is evaluated independently at each step. The number of steps before such a relocation takes place is finite however small $q_1$ may be as long as it remains strictly positive. The number of steps N to the first relocation has a geometric distribution with success probability equal to $q_1$ with mass distribution function

$$prob\left\{N=k\right\}=q_1{\left(1-q_1\right)}^k,k\ge 0.$$

(34)

Therefore, the expected value and variance of N, respectively, are

$$\mathrm{E}\left[N\right]=\frac{1}{q_1}\mathrm{and}Var\left[N\right]=\frac{1-q_1}{q_1^2}.$$

(35)

Even though the relocation modification clearly violates the inherent locality of the Walktrap family of heuristics, if properly calibrated, it happens infrequently enough so as not to severely degrade time performance. Moreover, in a distributed system, a simple move of the random walker to the appropriate graph segment suffices and its cost is certainly affordable. Moreover, since relocation is a rare event, the total number of relocations can be modeled by a Poisson distribution.

5. Analysis

5.1. Data

In order to experimentally evaluate the performance of Markov Walktrap, a Kronecker synthetic graph [59,60,61] has been created. Kronecker graphs are recursively constructed from an original generator graph with the following model

$$\begin{array}{cc}\hfill {\mathbf{A}}_0& =\mathbf{A},\hfill \\ \hfill {\mathbf{A}}_{k+1}& ={\mathbf{A}}_k\otimes \mathbf{A},n\ge 1,\hfill \end{array}$$

(36)

where ${\mathbf{A}}_0$ is the generator graph and ⊗ denotes the Kronecker tensor product.

The generator matrix was

$${\mathbf{A}}_0=\left[\begin{array}{ccccccc}1& 1& 1& 0& 0& 1& 1\\ 1& 1& 1& 1& 0& 1& 1\\ 1& 1& 1& 1& 1& 0& 0\\ 0& 1& 1& 1& 1& 1& 0\\ 0& 0& 1& 1& 1& 1& 1\\ 1& 1& 0& 1& 1& 1& 1\\ 1& 1& 0& 0& 1& 1& 1\end{array}\right],$$

(37)

which has $p=7$ vertices numbered from 0 to 6 and each $v_k$ is connected to $v_{k_1}$, $v_{k_2}$, and $v_{k_3}$, where

$$\begin{array}{c}{k}_1=\left(k-1\right)modp,\\ k_2=\left(k+1\right)modp,\\ k_3=\left(k+2\right)modp.\end{array}$$

(38)

The generator graph is shown in Figure 1.

The Kronecker model of Equation (36) has been executed six times aiming to obtain a large graph Y whose properties are summarized in

Table 2. Structural properties of graph Y.

Property	Value	Property	Value	Property	Value	Property	Value
Vertices	117,649	${\sigma }_0$	7.68 × 10${}^{-4}$	Triangles	37127	Squares	1981
Edges	5,315,625	${\sigma }_0^{\prime }$	$0.6631$	min cost	$5.1891$	max cost	$10.7232$
${\rho }_0$	$45.1821$	Diameter length	19	avg cost	$72.1149$	avg cost	$106.0012$
${\rho }_0^{\prime }$	$1.3263$	Diameter cost	$124.4021$	max cost	$143.2716$	max cost	$198.2221$

.

Definitions 5 and 6 outline four important structural graph metrics.

Definition 5.

The (log)density of a fuzzy graph is the ratio of the (logarithm of the) number of its edges to the (logarithm of the) number of its vertices.

$$\begin{array}{c}{\rho }_0\stackrel{▵}{=}\frac{\left|E\right|}{\left|V\right|}=\frac{{\left|V\right|}^{1+ϵ}}{{\left|V\right|}^{\phantom{1+ϵ}}}={\left|V\right|}^{ϵ},0\le ϵ\le 1,\hfill \\ {\rho }_0^{\prime }\stackrel{▵}{=}\frac{log\left|E\right|}{log\left|V\right|}=\frac{log{\left|V\right|}^{1+ϵ}}{log\left|V\right|}=1+ϵ.\hfill \end{array}$$

(39)

Definition 6.

The (log)completeness of a fuzzy graph is the ratio of the (logarithm of the) number of its edges to the (logarithm of the) number of the edges of the complete graph with the same number of vertices.

$$\begin{array}{c}{\sigma }_0\stackrel{▵}{=}\frac{\left|E\right|}{\left(\genfrac{}{}{0pt}{}{\left|V\right|}{2}\right)}=\frac{2\left|E\right|}{\left|V\right|\left(\left|V\right|-1\right)}\approx 2\frac{{\left|E\right|}^{\phantom{2}}}{{\left|V\right|}^2}=2{\left|V\right|}^{ϵ-1},\hfill \\ {\sigma }_0^{\prime }\stackrel{▵}{=}\frac{log\left|E\right|}{log\left(\genfrac{}{}{0pt}{}{\left|V\right|}{2}\right)}=\frac{log\left|E\right|}{log\left(\frac{\left|V\right|\left(\left|V\right|-1\right)}{2}\right)}=\frac{log\left|E\right|}{log\left|V\right|+log\left(\left|V\right|-1\right)-1}\approx \frac{log\left|E\right|}{2log\left|V\right|}=\frac{{\rho }_0^{\prime }}{2}=\frac{1+ϵ}{2}.\hfill \end{array}$$

(40)

Observe that Y is highly connected as it has a low diameter of a relatively low cost, high average degree, and a high number of triangles and squares.

5.2. Time and Memory Requirements

In

Table 3. Performance in terms of time (sec) and vertex visits.

Algorithm	Distribution	RW (s)	CB (s)	Total (s)	Visits
CW	Poisson	$128.1381$	$2000.9831$	$2129.1212$	471,858
CW	Binomial	$131.2182$	$2000.0304$	$2131.2486$	470,342
CW + I	Poisson	$144.4752$	$1911.9118$	$2056.3870$	477,423
CW + I	Binomial	$139.9916$	$1904.0216$	$2044.0132$	477,216
CW + R	Poisson	$157.0104$	$1840.0013$	$1997.0117$	476,997
CW + R	Binomial	$157.6633$	$1850.1003$	$2007.7636$	476,957
CW + IR	Poisson	$164.3779$	$2002.7745$	$2167.1524$	477,762
CW + IR	Binomial	$162.0222$	$1911.8664$	$2073.8886$	477,002
MW	Poisson	$157.5248$	$2104.9918$	$2262.5166$	475,308
MW	Binomial	$156.9924$	$2099.5256$	$2256.5180$	475,121
MW + I	Poisson	$165.3374$	$1908.6612$	$2073.9986$	478,313
MW + I	Binomial	$163.0015$	$1902.4319$	$2065.4334$	478,102
MW + R	Poisson	$173.9016$	$1850.0971$	$2023.9987$	477,916
MW + R	Binomial	$171.0017$	$1840.0054$	$2011.0071$	477,831
MW + IR	Poisson	$181.0013$	$2000.9917$	$2181.9930$	478,514
MW + IR	Binomial	$178.0017$	$2000.8585$	$2178.8602$	478,333
FW	Poisson	$191.9989$	$2425.1121$	$2617.1110$	515,444
FW	Binomial	$184.0451$	$2417.3376$	$2601.3827$	514,312
FN-G	Poisson	–	–	$9322.9514$	–
FN-G	Binomial	–	–	$9344.7778$	–

, the total execution time for Chebyshev Walktrap (CW), Markov Walktrap (MW), Fuzzy Walktrap (FW), and Fuzzy Newman–Girvan (FN-G) is shown. The last two algorithms are outlined in [10], whereas Fuzzy Markov was proposed in [15]. The effect of the escape mechanisms of weight inversion (I) and relocation (R) are also shown for Markov Walktrap and Chebyshev Walktrap. The Fuzzy Newman–Girvan is an exhaustive algorithm that will serve as baseline both for the requirements and the clustering quality. For the Walktrap algorithms, the time for the two phases, namely random walking (RW) and community building (CB), are recorded separately, while for the Fuzzy Newman–Girvan case only, the total time is recorded as there is only a single phase. In addition, the last column of Table 3 lists the number of the vertices visited by the random walker.

Fuzzy Newman–Girvan is considerably slower than any member of the Walktrap family of algorithms. This can be attributed to the exhaustive nature of Fuzzy Newman–Girvan as well as to the extensive use of locality by the Walktrap family. Moreover, the probabilistic constraints of Markov Walktrap and Fuzzy Walktrap resulted in the acceleration of both phases of the respective algorithms, with the second order constraints yielding the lowest times in each case. Concerning the escape strategy of the random walker, the relocation option resulted in a slower walking phase but in an accelerated community building phase, with that combination being more efficient than both weight inversion and the combination of the two escape strategies. Omitting an escape strategy is not advisable. Therefore, it is not recommended to activate both escape strategies at the same time. At any rate, the Chebyshev Walktrap with relocation (CW+R) had the best overall performance tagged along in a close manner by the Markov Walktrap with relocation (MW+R). The original Fuzzy Markov being the tardiest member of the family.

An explanation for the time achieved under the relocation strategy is that the teleportation of the random walker results in cache misses, which translates to expensive fetch cycles in the memory hierarchy system. This can be seen in the last two columns of Table 3, as there is not a clear correspondence between the number of total visits and the total walking phase time. When relocation is enabled, the mean visit time is clearly higher. At any rate, the number of visits is linear in the vertex set cardinality.

In addition, the selection of h did not appear to have a significant performance impact, although in most cases the random walker was slower when h was a Poisson random variable both in terms of time and in terms of total visits. This can be attributed to the large number of high cost edges which forced the walker to bounce more times inside a community before eventually moving to another. On the other hand, the symmetric form of the binomial distribution mass function resulted in a larger number of low cost edges, facilitating the movement of the random walker and making the communities easily separable compared to the Poisson case.

The memory requirements were monitored with the Ubuntu Watch administrative tool as presented in

Table 4. Performance in terms of memory (rounded in MBs).

Algorithm	Distribution	min	max	mean	std
CW	Poisson	4128	6742	4892	322
CW	Binomial	4128	6744	4880	324
CW + I	Poisson	4128	6739	4886	335
CW + I	Binomial	4128	6744	4881	338
CW + R	Poisson	4128	8002	5113	427
CW + R	Binomial	4128	8000	5121	422
CW + IR	Poisson	4128	8001	5208	345
CW + IR	Binomial	4128	8002	5214	339
MW	Poisson	4128	6744	4800	331
MW	Binomial	4128	6740	4881	325
MW + I	Poisson	4128	6739	4879	336
MW + I	Binomial	4128	6751	4883	335
MW + R	Poisson	4128	8004	5112	428
MW + R	Binomial	4128	8002	5108	428
MW + IR	Poisson	4128	8002	5200	343
MW + IR	Binomial	4128	8000	5204	341
FW	Poisson	4128	6750	4912	281
FW	Binomial	4128	6742	4910	280
FN-G	Poisson	8192	12,402	11,012	278
FN-G	Binomial	8192	12,464	11,121	280

. In contrast to other similar tools such as htop, Watch generates a text output which can be parsed and analyzed. It was periodically ran every 10 s through a bash script resulting in records of several thousand of entries each.

Fuzzy Newman–Girvan consumes more memory than any other algorithm presented in this article by far. However, it utlilizes the memory constantly and consistently, as denoted by the relatively low standard deviation. This is an important performance feature for operating systems process schedulers [46]. On the other hand, the Walktrap family exploits graph caching. This in turn translates to lower traffic between the disk and the memory, as Neo4j is not a memory-only database, as well as to fewer synchronization and serialization operations. When the relocation strategy is selected, then memory utilization has certain spikes, as it can be inferred from the increased maximum memory occupied and the increased standard deviation. This is a direct result of the random walker teleportation which temporarily annuls any scheduling optimization as well as any caching done at the software or hardware level.

5.3. Community Coherence

The following definition will facilitate further analysis of the experimental results.

Definition 7.

The (log)scree plot of a set S is the plot of the (logarithm of the) values of S versus their sorted frequency.

Since Y does not contain ground truth communities, the communities obtained by the Fuzzy Newman–Girvan will be used as a baseline reference since their sizes are closer to a power law distribution, which is an essential feature of large, scale-free graphs. The deviation $\xi$ of a set of numbers ${\left\{x\right\}}_{k=1}^n$ from a power law

$$f_k={\alpha }_0{k}^{-{\gamma }_0},1\le k\le n,{\alpha }_0>0,{\gamma }_0>0$$

(41)

is quantified by the formula [62,63]

$$\xi =\sqrt{\frac{1}{n}\sum _{k=1}^n{\left(log{f}_k-\left(log{\alpha }_0-{\gamma }_{0}logk\right)\right)}^2},$$

(42)

where parameters ${\alpha }_0$ and ${\gamma }_0$ can be estimated by, for instance, a least squares method [25]. Additionally, the estimated value of ${\alpha }_0$ serves as a quality indicator, as it should be as close to $[2,3]$ as possible.

The number of communities for each algorithm are shown in

Table 5. Number of communities.

	CW	CW + I	CW + R	CW + IR	MW	MW + I	MW + R	MW + IR	FW	FN-G
Poisson	13,751	137,58	13,332	13,443	13,761	13,789	13,456	13,804	14,127	12,816
Binomial	18,841	18,912	16,090	17,621	18,877	18,801	17,002	18,811	18,891	15,117

. Notice that this is not an absolute clustering quality metric, as typically a large number of coherent communities is preferable to a smaller number of sparse ones. Nonetheless, the introduction of the relocation strategy systematically pushes the number of communities towards the reference number, although more evidence is required for determining community coherence. This will be addressed by the two asymmetric indices of this section.

In order to evaluate the clustering quality, the Kullback–Leibler divergence between the sorted sizes of the communities generated by the Fuzzy Newman–Girvan and the sorted community sizes of the remaining algorithms was computed. Recall that for two discrete distributions $p_k$ and $q_k$ the Kullback–Leibler divergence is defined as

$$〈p\left|\right|q〉\stackrel{▵}{=}\sum _{k=1}^n{p}_{k}log\left(\frac{p_k}{q_k}\right)=\sum _{k=1}^n{p}_{k}log{p}_k-\sum _{k=1}^n{p}_{k}log{q}_k,$$

(43)

where k ranges over the union of discrete events. If $p_k$ and $q_k$ have no events in common, then the result is undefined. If for a single event $p_k=0$ or $q_k=0$, then the corresponding summand is zero.

Table 6. Kullback–Leibler divergence.

	CW	CW + I	CW + R	CW + IR	MW	MW + I	MW + R	MW + IR	FW
Poisson	$0.6409$	$0.6311$	$0.2766$	$0.4504$	$0.3826$	$0.3911$	$0.3281$	$0.4519$	$0.8012$
Binomial	$0.3885$	$0.3977$	$0.3519$	$0.3700$	$0.5804$	$0.5687$	$0.6140$	$0.5626$	$0.6748$

summarizes the divergence for the Poisson and the binomial cases.

Chebyshev Walktrap with relocation outperforms the remaining algorithms, as it has less divergence from the reference distribution.

A question at this point is whether a correspondence between the communities returned by each algorithm can be found. The asymmetric Tversky index between two sets T and V is defined as

$${\nu }_{T,V}\stackrel{▵}{=}\frac{\left|T\cap V\right|}{\left|T\cap V\right|+w_1\left|T\setminus V\right|+w_2\left|V\setminus T\right|},w_1,w_2>0,$$

(44)

and it quantifies the distance between the template set T and the variant set V. By the very definition of the index, the template set T and the variant set V are not interchangeable, namely ${\nu }_{T,V}\ne {\nu }_{V,T}$. This agrees with intuition, as it makes sense to ask how much the heuristic results differ from the ground truth community, whereas there is no point in asking the inverse question. On the contrary, with a symmetric distance metric such as, for instance, the Tanimoto similarity coefficient

$${\tau }_{T,V}\stackrel{▵}{=}\frac{\left|T\cap V\right|}{\left|T\cup V\right|}=\frac{\left|T\cap V\right|}{\left|T\right|+\left|V\right|-\left|T\cap V\right|},$$

(45)

no distinction can be made between the template and the variant, which can potentially lead to misleading results.

At this point, it should be highlighted that Fuzzy Newman Girvan was executed only once since it is a deterministic algorithm.

Returning to Label (44), the case $w_1+w_2=1$ is of particular interest in data mining, as it confines the coefficients on the plane which maximizes the minimum distance of T from V. Notice that algebraically this asymmetry stems from both the terms $\left|T\setminus V\right|$ and $\left|V\setminus T\right|$, which denote the number of elements of T not found in V and vice versa. Both terms signify in their own way how V is different from T. The former corresponds to the part of V which is missing from T, whereas the latter corresponds to any additions to V. As a rule, $\left|T\setminus V\right|$ is more important and, consequently, $w_1>w_2$. As there is no standard rule for selecting $w_1$ and $w_2$, the following two schemes have been used, a linear

$$w_1=\frac{s}{1+s}=\frac{1}{1+\frac{1}{s}},w_2=\frac{1}{1+s},1\le s\le 5$$

(46)

and an exponential

$$w_1=\frac{e^s}{1+e^s}=\frac{1}{1+e^{-s}},w_2=\frac{1}{1+e^s},0\le s\le 2.$$

(47)

Observe that in the first case $\frac{w_1}{w_2}=s$, while in the second $\frac{w_1}{w_2}=e^s$, which clearly represents a non-linear scaling of the first case. Furthermore, the second case is considerably biased in favor of $\left|T\setminus V\right|$.

Once for each possible pair of the m ground truth communities $T_i$$1\le i\le m$ and the n estimated ones $V_j$$1\le j\le n$ the $mn$ Tversky indices have been computed, the similarity score $J\left(s\right)$ for a given s is computed

$$J\left(s\right)\stackrel{▵}{=}\frac{1}{mn}\sum _{i=1}^m\sum _{j=1}^n{\nu }_{T_i,V_j}.$$

(48)

Summing over the range of s and taking the average, the mean similarity score $\overline{J}$ is obtained

$$\overline{J}\stackrel{▵}{=}\frac{1}{\left|s\right|}\sum _{s}J\left(s\right).$$

(49)

The overall similarity scores are shown in

Table 7. Tversky index.

Linear	CW	CW + I	CW + R	CW + IR	MW	MW + I	MW + R	MW + IR	FW
Poisson	$0.6199$	$0.6126$	$0.6616$	$0.6012$	$0.5881$	$0.5902$	$0.6211$	$0.5577$	$0.2742$
Binomial	$0.4323$	$0.4153$	$0.6059$	$0.5853$	$0.5648$	$0.5702$	$0.5814$	$0.5332$	$0.2131$
Exponential	CW	CW + I	CW + R	CW + IR	MW	MW + I	MW + R	MW + IR	FW
Poisson	$0.5505$	$0.5345$	$0.7503$	$0.5462$	$0.5271$	$0.5383$	$0.7354$	$0.4811$	$0.1703$
Binomial	$0.5230$	$0.4841$	$0.6992$	$0.4737$	$0.5102$	$0.4890$	$0.6772$	$0.4283$	$0.1523$

.

Again, Chebyshev Walktrap with relocation outperforms the remaining algorithms as it has the highest similarity with the reference communities. Note that the exponential weighting scheme sharpens the difference between the algorithms by raising the maximum scores and lowering the minimum ones.

For the experiments of the section, the termination criterion ${\tau }_0$ was chosen to be a user supplied number of iterations, namely $\left|V\right|log\left|V\right|$. This number of iterations is sufficiently large for generating communities in a reliable way. Moreover, each iteration is very quick, so the overall execution time was kept at an acceptable level despite the large number of iterations.

5.4. Relocations

Analysis is concluded with a summary regarding the relocations made by the Chebyshev Walktrap and the Markov Walktrap.

In

Table 8. Relocation summary.

	CW + R	MW + R
Number of relocations	18	14
First relocation step	101	44
min between relocations	17,812	32,991
max between relocations	27,099	64,818
mean between relocations	21,002	38,002

, certain statistics regarding the random walker relocations are shown. Specifically, the first line presents the total number of relocations, whereas the second line shows the number of steps that the random walker makes before being relocated for the first time. Similarly, the last three lines contain the minimum, maximum and average number of steps between two successive relocations, respectively.

6. Conclusions

The primary contribution of this article is the implementation over Neo4j of Chebyshev Walktrap, a community discovery algorithm designed for edge-fuzzy graphs, a class of fuzzy graphs used among others in [10], which is based on the Random Walker algorithmic principle. Additionally, Chebyshev Walktrap relies on the competitive factors of second order statistics though the Chebyshev inequality and on an optional relocation capability in order to bound unnecessarily costly walks and, thus, remaining inside a community and being trapped for too long within the boundaries of a community, respectively. The relocation aspect was also backported to the Markov Walktrap algorithm first proposed in [15]. The effect of relocation on the community coherence was evaluated based on the asymmetric Tversky index using the Fuzzy Newman–Girvan algorithm from [15] as baseline, while its effect on the output distribution was assessed with the asymmetric Kullback–Leibler divergence. The latter was also the basis for evaluating the distance between the community size distribution generated by Fuzzy Newman–Girvan and the one computed by the Makrov Walktrap and the Chebyshev Walktrap. In these cases, the introduction of asymmetry resulted in the clear distinction between the baseline data and their variants.

The test dataset was a large synthetic Kronecker graph whose edge fuzziness was controlled either by a binomial or by a Poisson distribution. In this dataset, our performance metrics showed that Chebyshev Walktrap yields more compact communities whose sizes are more clustered. Additionally, Markov Walktrap is, in many instances, slightly faster at the expense of a somewhat bigger memory footprint.

The experimental results of Section 5 hint at some future research directions. More sophisticated from a probabilistic viewpoint, community discovery algorithms should be able to exploit the asymmetry of the edge fuzziness distribution through higher order concentration inequalities such as the Talagrand inequality, provided their computation is efficient. Moreover, new metrics for community matching, perhaps utilizing functional or semantic information should be developed. Additionally, methodologies for reliably assessing community coherence based on higher order structural or functional interactions should be sought. Finally, more experiments in larger graphs should be conducted in order to determine any inherent scalability limitations.