Preferential placement for community structure formation

TL;DR

This paper introduces the principle of preferential placement to model realistic community structures in complex networks, supported by empirical analysis and theoretical insights.

Contribution

It proposes a novel principle for modeling community structures and provides both empirical and theoretical analysis of the resulting network properties.

Findings

Community structures can be effectively modeled using preferential placement.

The empirical analysis shows realistic clustering patterns.

Theoretical results support the validity of the proposed model.

Abstract

Various models have been recently proposed to reflect and predict different properties of complex networks. However, the community structure, which is one of the most important properties, is not well studied and modeled. In this paper, we suggest a principle called "preferential placement", which allows to model a realistic clustering structure. We provide an extensive empirical analysis of the obtained structure as well as some theoretical results.

Figures

Equations

$$\mathrm{E}{F}_n(s)=\frac{c(s-1)!\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)}{(2-c)\mathrm{\Gamma }\left(s+1+\frac{1}{1-c}\right)}\left(n+O\left(s^{\frac{1}{1-c}}\right)\right)\sim \frac{c\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)}{(2-c)}\cdot \frac{n}{s^{1+\frac{1}{1-c}}}.$$

$$\mathrm{E}{F}_n(s)=\frac{c(s-1)!\mathrm{\Gamma }(3-\alpha )}{(2-\alpha )\mathrm{\Gamma }(s+2-\alpha )}\left(n^{1-\alpha }+O\left(n^{\max \{0,1-2\alpha \}}{s}^{1-\alpha +ϵ}\right)\right)\sim \frac{c\mathrm{\Gamma }(3-\alpha )}{2-\alpha }\cdot \frac{n^{1-\alpha }}{s^{2-\alpha }}.$$

$$X_{\delta }(n)=\frac{1}{n}\sum _{i=1}^n\Vert {\mathbf{x}}_i{\Vert }^{\delta },\delta <\beta ,$$

$$\mathrm{E}{X}_{\delta }(n)\le \mathrm{E}{\xi }^{\delta }\sum _{t=2}^n\frac{1}{t}\sim \mathrm{E}{\xi }^{\delta }\log n.$$

$$\mathrm{E}{X}_2(n)=\mathrm{E}{\xi }^2\sum _{t=2}^n\frac{1}{t}\sim \mathrm{E}{\xi }^2\log n.$$

$$\mathrm{E}\Vert \mathbf{\xi }+\mathbf{\eta }{\Vert }^{\delta }\le \mathrm{E}\Vert \mathbf{\xi }{\Vert }^{\delta }+\mathrm{E}\Vert \mathbf{\eta }{\Vert }^{\delta }.$$

$$\begin{gathered} \mathrm{E}{X}_2(n+1)=\frac{n}{n+1}\mathrm{E}{X}_2(n)+\frac{1}{n+1}\cdot \frac{1}{n}\cdot \sum _{i=1}^n\mathrm{E}\Vert {\mathbf{x}}_i+\mathbf{\xi }{\Vert }^2 \\ =\frac{n}{n+1}\mathrm{E}{X}_2(n)+\frac{1}{n+1}\cdot \frac{1}{n}\cdot \sum _{i=1}^n\mathrm{E}(\Vert {\mathbf{x}}_i{\Vert }^2+\Vert \mathbf{\xi }{\Vert }^2+2{\mathbf{x}}_i\cdot \mathbf{\xi })=\mathrm{E}{X}_2(n)+\frac{1}{n+1}\mathrm{E}\Vert \mathbf{\xi }{\Vert }^2. \end{gathered}$$

$$\mathrm{E}{X}_2(n)=\mathrm{E}\Vert \mathbf{\xi }{\Vert }^2\cdot \sum _{t=2}^n\frac{1}{t}\sim \mathrm{E}{\xi }^2\cdot \log n.$$

$$\begin{gathered} \mathrm{E}{X}_{\delta }(n+1)=\frac{n}{n+1}\mathrm{E}{X}_{\delta }(n)+\frac{1}{n+1}\cdot \frac{1}{n}\cdot \sum _{i=1}^n\mathrm{E}\Vert {\mathbf{x}}_i+\mathbf{\xi }{\Vert }^{\delta }\le \frac{n}{n+1}\mathrm{E}{X}_{\delta }(n) \\ +\frac{1}{n+1}\cdot \frac{1}{n}\cdot \sum _{i=1}^n\mathrm{E}\Vert {\mathbf{x}}_i{\Vert }^{\delta }+\frac{1}{n+1}\cdot \frac{1}{n}\cdot \sum _{i=1}^n\mathrm{E}\Vert \mathbf{\xi }{\Vert }^{\delta }=\mathrm{E}{X}_{\delta }(n)+\frac{1}{n+1}\mathrm{E}\Vert \mathbf{\xi }{\Vert }^{\delta }, \end{gathered}$$

$$X_{max}(n)=\underset{1\le i\le n}{\max }\Vert {\mathbf{x}}_i\Vert .$$

$$\mathrm{P}\left(X_{max}(n)>\frac{n^{\frac{1}{\beta }}}{\omega (n)}\right)=1-o(1).$$

$$\begin{gathered} \mathrm{P}(X_{max}(n)\le r)\le \mathrm{P}({\xi }_t\le r\text{ for all 1≤t≤n})={\left(1-\frac{1}{(1+r{)}^{\beta }}\right)}^n \\ \sim e^{\frac{-n}{(1+r{)}^{\beta }}}=o(1). \end{gathered}$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}=I\big{[}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq\theta\big{]}\,,\,\,\,\,\,\,\theta>0\,;$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}=pI\big{[}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq\theta\big{]}\,,\,\,\,\,\,0<p<1,\,\theta>0\,;$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}=p_{0}+p_{1}I\big{[}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq\theta\big{]}\,,\,\,\,\,\,0<p_{0},p_{1},p_{0}+p_{1}<1,\,\theta>0\,;$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}\propto\frac{1}{\|\mathbf{x}_{i}-\mathbf{x}_{j}\|}\,;$$

$$P\big{(}(v_{i},v_{j})\in E\big{)}\propto e^{-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|/d}\,,\,\,\,\,\,d>0\,.$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}=F(\|\mathbf{x}_{i}-\mathbf{x}_{j}\|,d_{tree}(v_{i},v_{j}))\,,$$

$$\mathrm{P}\big{(}(v_{i},v_{j})\in E\big{)}=I\big{[}\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq\theta\text{ or }d_{tree}(v_{i},v_{j}))=1\big{]}\,.$$

$$\underset{n\to \infty }{\lim }\frac{H_n}{\log n}=e.$$

$$\underset{n\to \infty }{\lim }\frac{\mathbf{diam}(G_n)}{\log n}\le 2e.$$

$$\mathrm{E}(F_{t+1}(1)|S_t)=F_t(1)\left(1-\frac{1-p(t)}{t}\right)+p(t),$$

$$\mathrm{E}(F_{t+1}(s)|S_t)=F_t(s)\left(1-\frac{s(1-p(t))}{t}\right)+F_t(s-1)\frac{(s-1)(1-p(t))}{t},s>1.$$

$$\mathrm{E}{F}_n(s)=\frac{c(s-1)!\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)}{(2-c)\mathrm{\Gamma }\left(s+1+\frac{1}{1-c}\right)}\left(n+{\theta }_{n,s}\right).$$

$$\mathrm{E}{F}_n(1)=\frac{c}{2-c}\left(n+{\theta }_{n,1}\right).$$

$$\begin{gathered} \mathrm{E}(F_{t+1}(1))=\mathrm{E}{F}_t(1)\left(1-\frac{1-c}{t}\right)+c=\frac{c}{2-c}\left(t+{\theta }_{t,1}\right)\left(1-\frac{1-c}{t}\right)+c \\ =\frac{c}{2-c}\left(t+1+{\theta }_{t,1}\left(1-\frac{1-c}{t}\right)\right). \end{gathered}$$

$$C\left(1-\frac{1-c}{t}\right)\le C,$$

$$\begin{gathered} \mathrm{E}(F_{t+1}(s))=\mathrm{E}{F}_t(s)\left(1-\frac{s\left(1-c\right)}{t}\right)+\mathrm{E}{F}_t(s-1)\frac{(s-1)\left(1-c\right)}{t} \\ =\frac{c(s-1)!\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)\left(t+{\theta }_{t,s}\right)}{(2-c)\mathrm{\Gamma }\left(s+1+\frac{1}{1-c}\right)}\left(1-\frac{s(1-c)}{t}\right)+\frac{c(s-1)!\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)(1-c)(t+{\theta }_{t,s-1})}{(2-c)\mathrm{\Gamma }\left(s+\frac{1}{1-c}\right)t} \\ =\frac{c(s-1)!\mathrm{\Gamma }\left(2+\frac{1}{1-c}\right)}{(2-c)\mathrm{\Gamma }\left(s+1+\frac{1}{1-c}\right)}\left(t+1+{\theta }_{t,s}\left(1-\frac{s(1-c)}{t}\right)+{\theta }_{t,s-1}\frac{s(1-c)+1}{t}\right). \end{gathered}$$

$$(s-1{)}^{\frac{1}{1-c}}\frac{s(1-c)+1}{t}\le s^{\frac{1}{1-c}}\frac{s(1-c)}{t}.$$

$${\left(1-\frac{1}{s-1}\right)}^{\frac{1}{1-c}}\ge 1+\frac{1}{s(1-c)}.$$

Full text

1 Introduction

The evolution of complex networks attracted a lot of attention in recent years. Empirical studies of different real-world networks have shown that such structures have some typical properties: small diameter, power-law degree distribution, clustering structure, and others [10, 17, 46]. Therefore, numerous random graph models have been proposed to reflect and predict such quantitative and topological aspects of growing real-world networks [10, 12, 17, 49, 55].

The most extensively studied property of complex networks is their vertex degree distribution. For the majority of studied real-world networks, the portion of vertices of degree $d$ was observed to decrease as $d^{-\gamma}$, usually with $2<\gamma<3$ [6, 23, 47]. Such networks are often called scale-free. The most well-known approach to modeling scale-free networks is called preferential attachment. The main idea of this approach is that new vertices emerging in a graph connect to some already existing vertices chosen with probabilities proportional to their degrees. Preferential attachment is a natural process allowing to obtain a graph with a power-law degree distribution, and many random graph models are based on this idea, see, e.g., [11, 14, 28, 35, 60].

Another important characteristic of complex networks is their community (or clustering) structure, i.e., the presence of densely interconnected sets of vertices, which are usually called clusters or communities [24, 26]. Several empirical studies have shown that community structure of different real-world networks has some typical properties. In particular, it was observed that the cumulative community size distribution obeys a power law with some parameter $\lambda$.111Cumulative community size distribution is the function defined for each $x$ as the probability that a size of a randomly sampled cluster is at least $x$. For instance, [16] reports that $\lambda=1$ for some networks; [3] obtains either $\lambda=0.5$ or $\lambda=1$; [27] also observes a power law with $\lambda$ close to 0.5 in some range of cluster sizes; [50] studies the overlapping communities and shows that $\lambda$ is ranging between 1 and 1.6.

Community structure is an essential property of complex networks. For example, it highly affects the spreading of infectious diseases in social networks [29, 39], spread of viruses over computer networks [58], promotion of products via viral marketing [32], propagation of information [56], etc. Therefore, it is crucial to be able to model realistic community structures.

Unfortunately, the widely used preferential attachment model fails to provide desired clustering structure [12] and only a few random graph models are able to generate realistic clusters. Probably the most well-known model was suggested in [37] as a benchmark for comparing community detection algorithms. In this model, the distributions of both degrees and community sizes follow power laws with predetermined exponents. However, there are two drawbacks of this model. First, it does not explain the power-law distribution of community sizes, these sizes are just sampled from a power-law distribution at the beginning of the process. Second, a subgraph induced by each community is very similar to the configuration model [9], which does not allow to model, e.g., hierarchical community structure often observed in real-world networks [3, 16].

A weighted model which naturally generates communities was proposed in [36]. However, the community structure in this model is not analyzed in detail and only the local clustering coefficient is shown. From the figures presented in [36] it seems that the community size distribution does not have a heavy tail as it is observed in real-world complex networks.

Another promising model which naturally generates clusters is proposed in [5] and it is based on a so-called clustering attachment. However, again, the distribution of cluster sizes was not analyzed for this model. Although from the presented figures and simulations222http://rocs.hu-berlin.de/interactive/ca/index.html it seems that this distribution is not heavy tailed, it would be interesting to analyze it theoretically or empirically.

Let us also mention a paper [52] which analyzes a community graph, where vertices refer to communities and edges correspond to shared members between the communities. The authors show that the development of the community graph seems to be driven by preferential attachment. They also introduce a model for the dynamics of overlapping communities. Note that [52] only models the membership of vertices and does not model the underlying network.

Finally, note that a simplified way to analyze the clustering structure is to measure the clustering coefficient, i.e., the probability that two neighbors of a vertex are connected. Several models with high clustering coefficient were proposed in the literature. For example, Spatial Preferential Attachment model [2] was proven to have a power-law degree distribution, high clustering coefficient and other desirable properties [30, 53]. Similarly, hyperbolic random graphs [34] generate a power-law degree distribution, low diameter and a high clustering coefficient. Also, preferential attachment models were modified by introducing special steps of triangle formation [28]. However, a high clustering coefficient does not necessary lead to a realistic clustering structure. The presence of clusters and the distribution of their sizes were not analyzed in the models mentioned above.

In this paper, we propose a process which naturally generates clustering structures. Our approach is called preferential placement and it is based on the idea that vertices can be embedded in a multidimensional space of latent features. The vertices appear one by one and their positions are defined according to preferential placement: new vertices are more likely to fall into already dense regions. We present a detailed description of this process in Section 2. After $n$ steps we obtain a set of $n$ vertices placed in a multidimensional space. In Section 3 we empirically and theoretically analyze the obtained structure: in particular, we show that the communities are clearly visible and their sizes are distributed according to a power law. Note that after the placement of all vertices is defined, one can easily construct an underlying network, using, e.g., the threshold model [13, 41]. We discuss possible models and their properties in Section 4.

This paper is a journal version of [20]. It contains additional theoretical and empirical results for the vertex configuration obtained by preferential placement. We also discuss some new ideas on generating graphs from this configuration and analyze the properties of the obtained random graph models.

2 Preferential placement

2.1 Definition

In this section, we describe the proposed approach which we call preferential placement. We assume that all vertices are embedded in $\mathbb{R}^{d}$ for some $d\geq 1$. One can think that coordinates of this space correspond to latent features of vertices. Introducing latent features has recently become a popular approach both in predictive and generative models. These models are known by different names such as latent feature models [43, 44], matrix factorization models [4, 21, 42], spatial models [2, 7, 8], or geographical models [13, 41]. The basic idea behind all these models is that vertices having similar latent features are more likely to be connected by an edge.

Preferential placement is a procedure describing the embedding of vertices in the space $\mathbb{R}^{d}$. After that, given the coordinates of all vertices, one can construct a graph using one of many well-known approaches (see Section 4 for the discussion of possible variants).

Our model is parametrized by a distribution $\Xi$ taking nonnegative values. The proper choice of $\Xi$ is discussed further in this section.

We construct a random configuration of vertices (or points) $S_{n}=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$, where $\mathbf{x}_{i}=(x_{i}^{1},\ldots,x_{i}^{d})$ denotes the coordinates of the $i$-th vertex $v_{i}$. Let $S_{1}=\{\mathbf{x}_{1}\}$, $\mathbf{x}_{1}$ is the origin. Now assume that we have constructed $S_{t}$ for $t\geq 1$, then we obtain $S_{t+1}$ by adding a vertex $v_{t+1}$ with the coordinates $\mathbf{x}_{t+1}$ chosen in the following way:

• •

Choose a vertex $v_{i_{t+1}}$ from $v_{1},\ldots,v_{t}$ uniformly at random.

• •

Sample $\xi_{t+1}$ from the distribution $\Xi$.

• •

Sample a direction $\mathbf{e}_{t+1}$ from a uniform distribution on a multidimensional sphere $\|\mathbf{e}_{t+1}\|=1$, where by $\|\cdot\|$ we denote the Euclidean norm in $\mathbb{R}^{d}$.

• •

Set $\mathbf{x}_{t+1}=\mathbf{x}_{i_{t+1}}+\xi_{t+1}\cdot\mathbf{e}_{t+1}$.

We argue in this paper that in order to obtain a realistic clustering structure one should take $\Xi$ to be a heavy tailed distribution. In this case, according to the procedure described above, new vertices will usually appear in the dense regions, close to some previously added vertices; however, due to the heavy tail of $\Xi$, from time to time we get outliers, which originate new clusters.

We call the described above procedure “preferential placement” due to its analogy with preferential attachment. Assume that at some step of the algorithm we have several clusters, i.e., groups of vertices located close to each other, and a new vertex appears. Then the probability that this vertex will join a cluster $C$ is roughly proportional to its size, i.e., the number of vertices already belonging to this cluster. This is the basic intuition which we discuss further in this paper in more details.

2.2 Tree-based interpretation

In some cases, it is convenient to think of the definition given above in the following way. Let us first construct a uniform random recursive tree. Recall that a recursive tree of order $n$ is a rooted tree on $n$ vertices labeled $v_{1},\ldots,v_{n}$, with the property that for each $k$, $2\leq k\leq n$, the indices of the vertices on the unique path from the root $v_{1}$ to the vertex $v_{k}$ form an increasing sequence [57]. Uniform random recursive tree (or URRT) is a tree sampled from the set of recursive trees of order $n$ uniformly at random. An equivalent way to sample a URRT is to start from a graph $T_{1}$ consisting of a root vertex $v_{1}$ and at each step $t+1$ add a vertex $v_{t+1}$ connected to a vertex $v_{i_{t+1}}$ chosen from $v_{1},\ldots,v_{t}$ uniformly at random. After $n$ steps we get a URRT $T_{n}$. Random recursive trees were extensively studied in the literature, see, e.g., [18, 19, 45, 51].

To construct a random configuration of vertices $S_{n}=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$, we first construct a URRT $T_{n}$ using a recursive procedure described above (we will further call this tree genealogical). Then we label all edges of $T_{n}$: for each edge $(v_{t},v_{i_{t}})$ we set its label $\mathbf{w}_{t}$ to be a vector $\xi_{t}\mathbf{e}_{t}$, where $\xi_{t}$ is sampled from the distribution $\Xi$ and $\mathbf{e}_{t}$ is, as before, a random vector of length 1 in $\mathbb{R}^{d}$. Finally, the value $\mathbf{x}_{k}$ for each vertex $v_{k}$ is obtained by summing all labels on the unique path from $v_{1}$ to $v_{k}$, i.e., $\mathbf{x}_{k}=\mathbf{w}_{j_{1}}+\ldots+\mathbf{w}_{j_{l}}$ such that for all $1\leq t\leq l$ we have $i_{j_{t}}=j_{t-1}$, $j_{l}=k$, $j_{0}=1$. It is easy to see that this definition is equivalent to the one given in Section 2.1.

Note that each $\mathbf{x}_{k}=\mathbf{w}_{j_{1}}+\ldots+\mathbf{w}_{j_{l}}$ is obtained by a random walk with lengths of jumps distributed according to $\Xi$, which we further assume to be a power law. Such random walks are known in the literature as Lévy flights [15]. In other words, preferential placement is essentially a combination of URRT and Lévy flights.

3 Analysis of preferential placement

3.1 Experimental setup

In this section, we analyze structures obtained using the preferential placement procedure described above. We take $\Xi$ to be a slightly modified Pareto distribution with the density function $f_{\beta}(x)=\frac{\beta}{(x+1)^{\beta+1}},x\geq 0$ for fixed $\beta>0$. This distribution is supported on the interval $[0,+\infty)$ and is asymptotically equivalent to the Pareto distribution.

In all the experiments we take $d=2$ since the obtained structures are easy to visualize. However, we also tried other values of $d\geq 1$ and obtained results similar to shown, e.g., on Figures 5, 8 and 9. Also, if not specified otherwise, we generated structures with the number of points $n=10^{5}$. In Section 3.3.3, we analyze the effect of $d$, $n$ and $\beta$ on the distribution of cluster sizes.

3.2 Clustering structure depending on $\Xi$

First, let us visualize the structures obtained by our algorithm. We tried several values of $\beta$, $\beta\in\{0.5,1,1.5,2.5,4\}$. The results are presented on Figures 1 and 2. The value $\beta=0.5$ produces the heaviest tail, in this case the distribution $\Xi$ does not have a finite expectation. Although some clusters are clearly visible in this case (Figure 1(a)), they are located far apart from each other and there are too many single outliers laying far away from other points, which seems to be not very realistic. Note that for too large $\beta$, e.g., for $\beta=4$, the variance is too low and we obtain only one giant cluster with minor fluctuations, as presented on Figure 1(d). Further in this paper we discuss the case $\beta=1.5$ presented on Figure 2. In this case $\Xi$ has a finite expectation but an infinite variance.

Another interesting observation is a hierarchical clustering structure produced by our algorithm. To illustrate this, we take the figure obtained for $\beta=1.5$ and zoom it to see more details. Figure 2 shows that the largest cluster further consists of several sub-clusters.

In order to further analyze this self-similarity property, we empirically estimated a fractal dimension of the obtained structures [33]. A fractal dimension is an index characterizing how a fractal pattern changes with the scale at which it is measured. Unlike topological dimensions, the fractal dimension can take non-integer values. For example, a curve with a fractal dimension close to 1 behaves quite like an ordinary line, but a curve with a fractal dimension 1.9 winds convolutely through space nearly like a surface. We used a box-counting algorithm333https://en.wikipedia.org/wiki/Box_counting which computes the number of boxes $N(\varepsilon)$ with the edge size equal to $\varepsilon$ needed to cover the structure. Then the fractal dimension $d_{f}$ is defined using the following relation: $N(\varepsilon)\propto\varepsilon^{-d_{f}}$. The results are presented on Figure 3 (crosses represent the empirically observed pairs $(\varepsilon,N(\varepsilon))$, while lines correspond to the approximations $N(\varepsilon)\propto\varepsilon^{-d_{f}}$). One can see that the fractal dimension $d_{f}$ increases from $0.4$ to $1.4$ while $\beta$ increases from $0.5$ to $4$. Note that the fractal dimension of Lévy flights is also known to monotonically increase with the parameter of the power-law distribution [38].

3.3 The distribution of cluster sizes

In this section, we analyze the distribution of cluster sizes produced by preferential placement. We present both theoretical and empirical observations.

The main difficulty with the analysis of clustering structures is the fact that there are no standard definitions of clusters, both in graphs and metric spaces. For example, clusters are often defined as a result of some clustering algorithm.444Modularity, introduced in [48], can be used to define communities in graphs. However, this characteristic has certain drawbacks, as discussed in [25]. Moreover, modularity favors partitions with approximately equal communities, which contradicts the main idea of power-law distribution of community sizes. This causes a lot of difficulties for both theoretical and empirical analysis.

3.3.1 Theoretical analysis

First, let us discuss why we expect to observe a power-law distribution of cluster sizes in our model. As we discussed above, due to the absence of a rigorous definition of a cluster, further in this section we are only able to give some insights to the fact that the proposed algorithm is expected to generate power-law distribution of cluster sizes. Namely, we make some strong assumptions and then rigorously prove that the distribution follows a power law.

Let $F_{t}(s)$ denote the number of clusters of size $s$ at step $t$. In order to analyze $F_{t}(s)$ we consider its dynamics inductively. Assume that after a step $t$ we obtain some clustering structure. At step $t+1$ we add a vertex $v_{t+1}$ and choose its “parent” $v_{i_{t+1}}$ from $v_{1},\ldots,v_{t}$ uniformly at random. Clearly, the probability to choose a parent from some cluster $C$ with $|C|=s$ is equal to $\frac{s}{t}$. In this case, we call $C$ a parent cluster for $v_{t+1}$. Now let us make the following assumptions:

All clusters can only grow, they cannot merge or split. 2. 2.

At step $t+1$ a new cluster appears with probability $p(t)=\frac{c}{t^{\alpha}}$, $c>0$, $0\leq\alpha\leq 1$.555For simplicity, we choose $c$ such that $p(t)<1$ for all $t\geq 2$. Although our proof can be extended to cases where $p(t)$ behaves as $\frac{c}{t^{\alpha}}$ only asymptotically (which would allow to take larger $c$), we do not formally consider this case. 3. 3.

Given that a vertex $t+1$ does not create a new cluster, the probability to join a cluster $C$ with $|C|=s$ is equal to $\frac{s}{t}$.

These assumptions are quite strong and even not very realistic. For instance, it seems reasonable that two clusters can merge if many vertices appear somewhere between them. In Section 3.3.2 we discuss how the violation of the first assumption by the largest clusters affects the observed distribution of cluster sizes. Regarding the second assumption, $p(t)$ can possibly depend on the current configuration $S_{t}$. However, these assumptions allow us to analyze the behavior of $F_{t}(s)$ formally. Namely, we prove the following theorem.

Theorem 1

Under the assumptions described above the following holds.

If $\alpha=0$ and $0<c<1$, then

[TABLE] 2. 2.

If $0<\alpha\leq 1$, then for any $\epsilon>0$

[TABLE]

To sum up, if the probability $p(n)$ of creating a new cluster is of order $\frac{1}{n^{\alpha}}$ for $\alpha>0$, then the distribution of cluster sizes follows a power law with parameter $2-\alpha$ growing with $p(n)$ from 1 to 2; if $p(n)=c$, $0<c<1$, then the parameter grows with $c$ from 2 to infinity. Note that for power-law distributions, the parameter of the cumulative distribution $\lambda$ is one less than the parameter of the corresponding probability mass function, i.e., the exponent in the theorem discussed above is $\lambda+1$. The proof of Theorem 1 is technical and we place it in the Appendix.

Let us also explain why we do not consider $p(n)$ decreasing faster than $\frac{c}{n}$. It is natural to assume that a new cluster appears if a new vertex chooses a parent vertex near the border of some cluster and then $\xi_{t+1}$ and $\mathbf{e}_{t+1}$ are chosen such that $\mathbf{x}_{t+1}=\mathbf{x}_{i_{t+1}}+\xi_{t+1}\cdot\mathbf{e}_{t+1}$ falls quite away from the parent cluster. This probability is roughly proportional to the number of vertices located near the borders of the clusters. For each cluster, at least one vertex has to be near its border, so we get the lower bound $\frac{c}{n}$.

Finally, let us mention that in practice the probability $p(n)$ of creating a new cluster can depend not only on $\Xi$, but also on the definition of clusters. Further in this section we demonstrate that parameters of a clustering algorithm can affect the parameter of the obtained power law.

3.3.2 Empirical analysis.

As we already mentioned, there is no standard definition of a clustering structure. In many cases, clusters and communities are defined just as a result of some clustering algorithm. Therefore, we first analyze the performance of several clustering algorithms, then choose the most appropriate one and analyze clusters it produces.

We compare the following algorithms: k-means [40], EM (expectation maximization), and DBSCAN (density-based spatial clustering of applications with noise) [22]. For k-means and EM one has to specify the number of clusters. We tried several values of $k$, $k\in\{10,50,100,500,1000\}$, but both algorithms turned out to be not suitable for our problem. As expected, in all cases they unnaturally split the largest cluster into several small ones (see Figures 4(a) and 4(b)).

On the contrary, DBSCAN produces more realistic results. It requires two parameters: radius of neighborhood $\varepsilon$ and the minimum number of neighbors required to form a dense region $l$. We consider $l\in\{1,2,3\}$ and $\varepsilon$ is chosen in such a way that if we connect all pairs of vertices $i,j$ such that $\|i-j\|<\varepsilon$, then we get $Ln$ edges, $L\in\{5,25,125\}$, where $n$ is the number of vertices. For all parameters we get reasonable clustering structures. The result for $L=125,l=1$ is presented on Figure 4(c). For these parameters we also analyze the distribution of cluster sizes (see Figure 5(a)). Note that for not too large values of $s$ ($s<300$) the cumulative distribution follows a power law with parameter $\lambda\approx 0.95$. In Theorem 1 this value corresponds to the case $\alpha=0.05$, i.e., $p(n)\propto n^{-0.05}$. Based on this, we expect the number of clusters to grow as $n^{0.95}$, i.e., close to linearly. In Figure 6 we plot the empirical number of clusters and fit it by $n^{0.95}$.

Now, as we promised above, we show that $\lambda$ can depend on the clustering algorithm. Figure 5(b) shows the cumulative cluster size distribution for DBSCAN with $L=5,l=1$. We obtain $\lambda=1.44$, so it is larger in this case. Intuitively, the reason is that $p(n)$ is larger for $L=5$ than for $L=125$. Smaller values of $L$ correspond to smaller $\varepsilon$, which means that it is harder for a new vertex to join some existing cluster, which makes $p(n)$ larger.

Finally, let us also discuss the bend in the distribution observed on Figure 5. Namely, one can clearly see that several largest cluster sizes do not fit the desired power-law curve. This phenomenon is especially prominent for the largest cluster. The possible reason of this bend is the fact that in reality our first assumption in Section 3.3.1 is violated. Namely, we assumed that clusters cannot merge. However, Figure 4(c) clearly shows that the largest cluster further consists of several smaller ones (which is a desired hierarchical community structure). At previous steps of the construction process, these smaller clusters were separated, but then they merged. This fact caused the largest cluster to be larger than predicted, which we can see on Figure 5. Note that such phenomenon is expected only for a few largest clusters.

3.3.3 Influence of parameters

In this section, we discuss how the distribution of cluster sizes is affected by the parameters of the model. In particular, we show that our conclusions hold for various values of the number of vertices $n$, the parameter of Pareto distribution $\beta$, and the dimension $d$.

As a starting point, we take already discussed parameters $n=10^{5}$, $d=2$, $\beta=1.5$ and independently vary each parameter. The obtained results are presented in Figure 7.666According to the discussion above, for all curves, we removed the last point, which is an outlier corresponding to the largest cluster. We noticed that the distribution of cluster sizes is stable and is not affected much by the parameters of the model. The most different curve corresponds to the simple one-dimensional case $d=1$.

3.4 Scatter of configurations

In this section, we analyze the scatter of vertices in $\mathbb{R}^{d}$. This analysis is motivated by two factors. First, the scatter of points in $\mathbb{R}^{d}$ can be considered as an analogue of the diameter of graphs, which is an important and widely studied characteristic. Second, preferential placement is a simple and elegant process and its theoretical properties are of interest by themselves. Here we were partially motivated by the study of Lévy flights [15]: it is known that for Lévy flights we have $\langle|x|^{\delta}\rangle\propto t^{\delta/\beta}{\text{ if }}\delta\leq\beta$. Below we prove a similar result for preferential placement.

Recall that we assume the distribution $\Xi$ to have a density function $f_{\beta}(x)=\frac{\beta}{(x+1)^{\beta+1}}$, ${x\geq 0}$, ${\beta>0}$. Let us first consider the variable

[TABLE]

which is the average moment of order $\delta$ for the distance from the origin to a vertex. Recall that the most interesting case is $1<\beta<2$, for such $\beta$ the expected average squared distance $\mathrm{E}X_{2}(n)$ diverges, therefore we analyze $\mathrm{E}X_{\delta}(n)$ for $1\leq\delta\leq 2$. The following theorem holds.

Theorem 2

Let $\xi$ be a random variable sampled from the distribution $\Xi$. For any $\delta$, $1\leq\delta<\min\{\beta,2\}$:

[TABLE]

If $\beta>2$, then

[TABLE]

Note that in contrast to Lévy flights, preferential placement is much more concentrated near the origin (on average). Namely, $X_{\delta}(n)$ grows at most logarithmically with $n$ (instead of $n^{\delta/\beta}$). The reason is that preferential placement makes jumps from random vertices, while Lévy flight always chooses the latest one.

In the proof of Theorem 2 we will use the following lemma.

Lemma 1

Let $\boldsymbol{\xi}$ and $\boldsymbol{\eta}$ be two independent random vectors. Assume that the distribution of $\boldsymbol{\xi}$ is symmetrical, i.e., the distributions of $\boldsymbol{\xi}$ and $-\boldsymbol{\xi}$ are identical. Then for any $1\leq\delta\leq 2$:

[TABLE]

We place the proof of this lemma to Appendix. Now we are ready to prove Theorem 2.

*Proof. * Let us denote by $\boldsymbol{\xi}$ a vector of length $\xi$ (sampled from the distribution $\Xi$) pointing to a random direction, i.e., $\boldsymbol{\xi}=\xi\mathbf{e}$.

For $\beta>2$ the following recurrent formula holds:

[TABLE]

Using the above recurrent formula and the equality $\mathrm{E}X_{2}(1)=0$, we prove the second part of the theorem:

[TABLE]

To prove the first part we use Lemma 1:

[TABLE]

so we obtain the recurrent formula for $\mathrm{E}X_{\delta}(n)$ from which the theorem follows. $\Box$

In addition to $X_{\delta}(n)$, we analyze the following variable:

[TABLE]

Note that both $X_{\delta}(n)$ and $X_{max}(n)$ are interesting to analyze. For example, making parallels with graph theory, $X_{max}(n)$ is an analogue of the diameter, while $X_{\delta}(n)$ is similar to the average shortest path length [10, 46] which is often studied as a more stable analogue of the diameter.

For $X_{max}(n)$ the following simple lemma holds.

Lemma 2

Let $\omega(n)$ be any function tending to infinity as $n$ tends to infinity. Then

[TABLE]

*Proof. *

Let $r=\frac{n^{\frac{1}{\beta}}}{\omega(n)}.$ If for at least one step $t$ we have $\xi_{t}>2r$, then $X_{max}(n)>r$, therefore

[TABLE]

$\Box$

As a result, we obtain that “on average” preferential placement is well concentrated near the origin ($X_{\delta}(n)$ grows at most logarithmically). However, according to Lemma 2, the diameter of the obtained configuration grows at least as $n^{1/\beta}$, i.e., it is highly affected by outliers.

4 Graph models

4.1 Possible definitions

In this section, we discuss how a graph can be constructed based on the vertex embedding produced by the preferential placement procedure.

Spatial distance.

The basic idea behind many known spatial models is that we want to increase the probability of connecting two vertices if they have similar latent features. Various methods can be found in the literature, which are usually combined with some other ideas like introducing weights of vertices or taking into account degrees of vertices (see, e.g., a survey of spatial models in [8]). We now briefly describe some possible approaches:

• •

threshold model [13, 41]:

[TABLE]

• •

$p$-threshold model:

[TABLE]

• •

$p$-threshold model with random edges (as in spatial small-world models [8]):

[TABLE]

• •

inverted distance model:

[TABLE]

• •

Waxman model [59]:

[TABLE]

Here we denote by $E$ the set of edges and by $I[\cdot]$ the indicator function. We assume that all edges are mutually independent, hence to describe a random graph it is enough to define the probability of each edge. In Section 4.2 we present some empirical analysis for the basic threshold model as well as for its modification based also on a genealogical distance which we define in the next paragraph.

Genealogical distance.

It is also reasonable to take into account the genealogical tree used to create the configuration of vertices (see Section 2.2 for the description of this tree). For example, it seems very natural to connect each vertex to its parent in the genealogical tree. This serves several purposes: (i) the underlying genealogical process is reflected in the obtained graph, (ii) the graph is guaranteed to be connected, (iii) the graph has small diameter, as we show further in this section.

More generally, let $d_{tree}(v_{i},v_{j})$ be the length of the unique path between $v_{i}$ and $v_{j}$ in the genealogical tree. Then let

[TABLE]

where $F(x,y):\mathbb{R}_{+}\times\mathbb{N}\to[0,1]$ can be any function non-increasing with $x$ for fixed $y$ and with $y$ for fixed $x$.

In Section 4.2, in addition to the threshold model (1), we also analyze the following genealogical threshold model, which is a particular case of (2):

[TABLE]

4.2 Case study

As we promised above, in this section we analyze two models: the threshold model and the genealogical threshold model. In both cases, we choose $\theta$ such that we have $5n$ edges with $\|\mathbf{x}_{i}-\mathbf{x}_{j}\|\leq\theta$ in our graph. As before, we take $\Xi$ to be a distribution with the density function $f_{\beta}(x)=\frac{\beta}{(x+1)^{\beta+1}},x\geq 0$ for $\beta=1.5$.

Visualization.

Figure 8 visualizes the obtained graphs with $n=10^{4}$ (note that other empirical results are obtained for $n=10^{5}$, as usual). In both cases, one can clearly see communities of various sizes. One difference between two graphs is that in the threshold model the graph is not connected, as expected: 16% of vertices are isolated and the giant connected component consists of 54% of vertices.

Degree distribution.

Let us empirically analyze the degree distribution for both models (1) and (3). The cumulative degree distributions are presented on Figure 9. Observe that these distributions do not follow a power law. However, they are very similar to degree distributions obtained in many real-world networks (numerous examples can be found in [1]).

Distribution of community sizes.

Let us also empirically analyze the distribution of community sizes, which is the main focus of the current research. To do this, we partition the vertices using a well-known label propagation community detection algorithm (LPA) first proposed in [54]. The cumulative community size distributions are presented on Figure 10. Note that both distributions can be approximated by a power law. However, such approximation gives a better fit for the genealogical threshold model.

Diameter.

Now let us show that graphs obtained according to the genealogical threshold model (3) have small (at most logarithmic) diameter. In order to do this, we apply the following theorem [51].

Theorem 3

Let $H_{n}$ be the height of a uniform random recursive tree $T_{n}$. Then, with probability one,

[TABLE]

It remains to note that our graph contains the underlying genealogical tree by the definition, and this tree is URRT, so we obtain the following result.

Theorem 4

For a graph $G_{n}$ constructed according to (3) we have, with probability one,

[TABLE]

5 Conclusion and future work

In this paper, we introduced a principle called preferential placement. Our method is designed to model a realistic clustering structure in $\mathbb{R}^{d}$. The algorithm is parametrized only by a distribution $\Xi$, and if $\Xi$ is a Pareto distribution, which is a natural choice, then we essentially have only one parameter — the exponent $\beta$. The proposed algorithm naturally models clusters and the distribution of cluster sizes follows a power law, which is a desirable property. Although preferential placement only generates the coordinates of vertices, which is interesting in its own right, one can easily construct a graph based on the obtained structure using one of the methods discussed in this paper. In particular, we empirically analyzed the threshold model and its genealogical-based variant and obtained graphs with power-law distributions of community sizes and with degree distributions similar to ones observed in many real-worlds networks.

In this paper, we made only a first step to understanding the cluster formation in complex structures and there are many directions for future research. First of all, more formal analysis of the distribution of cluster sizes would be useful. As we discussed, the main problem here is the lack of any suitable formal definition of clusters. However, one can try, e.g., to analyze clusters produced by one of well-known clustering algorithms. Second, the radius of obtained configuration can be analyzed deeper: e.g., we are going to study the lower bound for $X_{\delta}(n)$ and the upper bound for $X_{max}(n)$ (see Section 3.4). Finally, there are many open questions in the analysis of the obtained graph structures: comparison of different models, theoretical analysis of the degree distribution, analysis of the local clustering coefficient, and so on.

Acknowledgements

This study was funded by RFBR according to the research project 18-31-00207 and by Russian President grant MK-527.2017.1.

Appendix

Proof of Theorem 1

First, recall the process of cluster formation:

• •

At the beginning of the process we have one vertex which forms one cluster;

• •

At step $t+1$ with probability $p(t)$ a new cluster consisting of $v_{t+1}$ is created;

• •

With probability $1-p(t)$ the new vertex $v_{t+1}$ joins some already existing cluster $C$ with probability proportional to $|C|$.

So, we can write the following equations:

[TABLE]

[TABLE]

Now we can take expectations of the both sides of the above equations and analyze the behavior of $\mathrm{E}F_{t}(s)$ inductively. Recall that the following form of the probability $p(t)$ is assumed: $p(t)=\frac{c}{t^{\alpha}}$, where $c>0$, $1\leq\alpha\leq 1$.

Proof for $\alpha=0.$

Consider the case $\alpha=0$, i.e., $p(t)=c$, $0<c<1$. Let us prove that in this case

[TABLE]

where $|\theta_{n,s}|\leq C\,s^{\frac{1}{1-c}}$ for some constant $C>0$.

We prove this result by induction on $s$ and for each $s$ the proof is by induction on $n$. Note that for $n=1$ Equation (6) holds for all $s$. Indeed, we can take $\theta_{1,s}=-1$ for $s>1$ and take $\theta_{1,1}$ satisfying $\frac{c}{2-c}\left(1+\theta_{1,1}\right)=1$.

Consider now the case $s=1$. We want to prove that

[TABLE]

For the inductive step we use Equation (4) and get

[TABLE]

Since

[TABLE]

this finishes the proof for $\alpha=0$ and $s=1$.

For $s>1$ we use Equation (5) and get

[TABLE]

To finish the proof we need to show that

[TABLE]

It is easy to show that the above inequality holds. Indeed, (7) can be rewritten as

[TABLE]

which holds as $\left(1-\frac{1}{s-1}\right)^{\frac{1}{1-c}}\geq 1-\frac{1}{(s-1)(1-c)}>1-\frac{1}{s(1-c)}$.

Proof for $0<\alpha\leq 1$.

Now we consider the case $p(t)=ct^{-\alpha}$ for $0<\alpha\leq 1$. Let us prove that in this case

[TABLE]

where $|\theta_{n,s}|\leq Cn^{\max\{0,1-2\alpha\}}s^{1-\alpha+\epsilon}$ for some constant $C>0$ and for any $\epsilon>0$.

The proof is similar to the case $\alpha=0$. Again, for $n=1$ the theorem holds, as we can take $\theta_{1,s}=-1$ for $s>1$ and take $\theta_{1,1}$ satisfying $\frac{c}{2-\alpha}\left(1+\theta_{1,1}\right)=1$.

Consider the case $s=1$. We want to prove that

[TABLE]

Inductive step in this case becomes

[TABLE]

The last equation holds since $(t+1)^{1-\alpha}=t^{1-\alpha}+(1-\alpha)t^{-\alpha}+O\left(t^{-\alpha-1}\right)$. In order to finish the proof for the case $s=1$ it is sufficient to show that

[TABLE]

which holds for sufficiently large $C$, since

[TABLE]

for all $\alpha$, $0<\alpha\leq 1$, and for sufficiently large $C$.

For $s>1$ we have:

[TABLE]

Recall the condition for $\theta_{n,s}$: $|\theta_{n,s}|\leq Cn^{\max\{0,1-2\alpha\}}s^{1-\alpha+\epsilon}$. Therefore, in order to finish the proof, we need to show that for some $C$ the sum of the error terms above can be bounded by $C(t+1)^{\max\{0,1-2\alpha\}}s^{1-\alpha+\epsilon}$, i.e.,

[TABLE]

Since $(t+1)^{\max\{0,1-2\alpha\}}\geq t^{\max\{0,1-2\alpha\}}$, it is sufficient to prove that

[TABLE]

which is equivalent to

[TABLE]

By using the inequality $(s-1)^{1-\alpha+\epsilon}\leq s^{1-\alpha+\epsilon}-(1-\alpha+\epsilon)s^{-\alpha+\epsilon}$, we get

[TABLE]

Note that the inequality

[TABLE]

holds for sufficiently large $C$, which completes the proof of the theorem.

Proof of Lemma 1

Let us first prove that for any vectors $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$ and for any $\delta$, $1\leq\delta\leq 2$, we have:

[TABLE]

If $\|\mathbf{x}\|=0$ and $\|\mathbf{y}\|=0$, then this inequality holds. Otherwise, without loss of generality, we can assume that $\|\mathbf{x}\|^{2}+\|\mathbf{y}\|^{2}=1$ (if $\|\mathbf{x}\|^{2}+\|\mathbf{y}\|^{2}=a>0$, we divide $\mathbf{x}$ and $\mathbf{y}$ by $\sqrt{a}$).

In order to prove (8), it is sufficient to show that $\frac{1}{2}\left(\|\mathbf{x}+\mathbf{y}\|^{\delta}+\|\mathbf{x}-\mathbf{y}\|^{\delta}\right)\leq 1$ and $1\leq\|\mathbf{x}\|^{\delta}+\|\mathbf{y}\|^{\delta}$. The second inequality is trivial:

[TABLE]

since $\|\mathbf{x}\|\leq 1,\|\mathbf{y}\|\leq 1$, and $\delta\leq 2$.

It remains to show that $\frac{1}{2}\left(\|\mathbf{x}+\mathbf{y}\|^{\delta}+\|\mathbf{x}-\mathbf{y}\|^{\delta}\right)\leq 1$. First let us note that

[TABLE]

Now we introduce the following variables

[TABLE]

and it remains to show that $2^{\delta/2-1}\left(u^{\delta/2}+v^{\delta/2}\right)\leq 1$, which is true since the maximum of $u^{\delta/2}+v^{\delta/2}=u^{\delta/2}+(1-u)^{\delta/2}$ for $0\leq u\leq 1$ is equal to $2^{-\delta/2+1}$ at $u=\frac{1}{2}$.

Now we note that the distributions of $\boldsymbol{\xi}+\boldsymbol{\eta}$ and $\boldsymbol{\xi}-\boldsymbol{\eta}$ are identical, so

[TABLE]

and now the lemma follows from Equation (8).