Factorization threshold models for scale-free networks generation

Many real networks such as the World Wide Web, financial, biological, citation and social networks have a power-law degree distribution. Networks with this feature are also called scale-free. Several models for producing scale-free networks have been obtained by now and most of them are based on the preferential attachment approach. We will offer the model with another scale-free property explanation. The main idea is to approximate the network's adjacency matrix by multiplication of the matrices $V$ and $V^T$, where $V$ is the matrix of vertices' latent features. This approach is called matrix factorization and is successfully used in the link prediction problem. To create a generative model of scale-free networks we will sample latent features $V$ from some probabilistic distribution and try to generate a network's adjacency matrix. Entries in the generated matrix are dot products of latent features which are real numbers. In order to create an adjacency matrix, we approximate entries with the Boolean domain $\{0, 1\}$. We have incorporated the threshold parameter $\theta$ into the model for discretization of a dot product. Actually, we have been influenced by the geographical threshold models which were recently proven to have good results in a scale-free networks generation. The overview of our results is the following. First, we will describe our model formally. Second, we will tune the threshold $\theta$ in order to generate sparse growing networks. Finally, we will show that our model produces scale-free networks with the fixed power-law exponent which equals two. In order to generate oriented networks with tunable power-law exponents and to obtain other model properties, we will offer different modifications of our model. Some of our results will be demonstrated using computer simulation.

adjacency matrix by a product of matrices V and V T , where V is the matrix of nodes' latent features vectors. To create a generative model of scale-free networks, we sample latent features V from some probabilistic distribution and try to generate a network adjacency matrix. Two nodes are connected by an edge if the dot product of their latent features exceeds some threshold. This threshold condition is influenced by the geographical threshold models that are applied to scale-free network generation [7]. Because of the methods used (adjacency matrix factorization and threshold condition), we call our model the factorization threshold model.
A network produced in such a way is scale-free and follows power-law degree distribution with an exponent of 2, which differs from the results for basic preferential attachment models [8][9][10] where the exponent equals 3. We also suggest an extension of our model that allows us to generate directed networks with a tunable power-law exponent.
This paper is organized as follows. "Related work" section provides information about related works that inspired us. The formal description of our model in the case of an undirected fixed size network is presented in "Model description" section, which is followed by a discussion of how to generate growing networks. In "Generating sparse networks" section, the problem of making resulting networks sparse is considered. "Degree distribution" section shows that our model indeed produces scale-free networks. Extensions of our model, which allows to generate directed networks with a tunable power-law exponents and some other interesting properties, will be discussed in "Model modifications" section. "Conclusion" section concludes the paper.

Related work
In this section, we consider related works that encouraged us to create a new model for complex networks generation.

Matrix factorization
Matrix factorization is a group of algorithms where a given matrix R is factorized into two smaller matrices Q and P such that: R ≈ Q T P [11].
There is a popular approach in recommendation systems which is based on matrix factorization [12]. Assume that users express their preferences by rating some items, this can be viewed as an approximate representation of their interests. Combining known ratings, we get partially filled matrix R, the idea is to approximate unknown ratings using matrix factorization R ≈ Q T P. A geometrical interpretation is the following. The rows of matrices Q and P can be seen as latent features vectors q i and p u of items and users, respectively. The dot product ( q i , p u ) captures an interaction between an user u and an item i, and it should approximate the rating of the item i by the user u: R ui ≈ (� q i , � p u ). Mapping of each user and item to latent features is considered as an optimization problem of minimizing distance between R and Q T P that is usually solved using stochastic gradient descent (SGD) or alternating least squares (ALS) methods.
Furthermore, matrix factorization was suggested to be used for link prediction in networks [6]. Link prediction refers to the problem of finding missing or hidden links which probably exist in a network [13]. In [6] it is solved via matrix factorization: a network adjacency matrix A is approximated by a product of the matrices V and V T , where V is the matrix of nodes' latent features.

Geographical threshold models
Geographical threshold models were recently proven to have good results in scale-free networks generation [7]. We are going to briefly summarize one variation of these models [14].
Suppose the number of nodes to be fixed. Each node carries a randomly and independently distributed weight variable w i ∈ R. Also, the nodes are uniformly and independently distributed with specified density in a R d . A pair of nodes with weights w, w ′ and Euclidean distance r are connected if and only if: where θ is the model threshold parameter and h(r) is the distance function that is assumed to decrease in r. For example, we can take h(r) = r −β , where β > 0.
First, exponential distribution of weights with the inverse scale parameter has been studied. This distribution of weights leads to scale-free networks with a power-law exponent of 2: P(k) ∝ k −2 . It is interesting that the exponent of a power law does not depend on the , d and β in this case. Second, Pareto weight distribution with scale parameter w 0 and shape parameter a has been considered. In this case, a tunable power-law degree distribution has been achieved: P(k) ∝ k −1− aβ d . There are other variations of this approach: uniform distribution of coordinates in the d-dimensional unit cube [15], lattice-based models [16,17] and even networks embedded in fractal space [18].

Model description
We studied theoretically matrix factorization by turning it from a trainable supervised model into a generative probabilistic model. When matrix factorization is used in machine learning, the adjacency matrix A is given and the goal is to train the model by tuning the matrix of latent features V in such way that A ≈ V T V . In our model, we make the reverse: latent features V are sampled from some probabilistic distribution and we generate a network adjacency matrix A based on V T V . Formally our model is described in the following way: • Network has n nodes and each node is associated with a d-dimensional latent features vector v i . • Each latent features vector v i is a product of weight w i and direction x i .
• Directions x i are i.i.d. random vectors uniformly distributed over the surface of (d − 1)-sphere. • Weights are i.i.d. random variables distributed according to Pareto distribution with the following density function f(w): • Edges between nodes i and j appear if a dot product of their latent features vectors ( v i , v j ) exceeds a threshold parameter θ.
Therefore, we take into consideration both node's importance w i and its location x i on the surface of a (d − 1)-sphere (that can be interpreted as the earth in the case of � x i ∈ S 2 ⊂ R 3 ). Thus, inspired by the matrix factorization approach we achieved the following model behavior: the edges in our model are assumed to be formed when a pair of nodes is spatially close and/or has large weights. Actually, compared with the geographical threshold models, we use dot product to measure proximity of nodes instead of Euclidean distance.
We have defined our model for fixed size networks, but in principle, our model can be generalized for the case of growing networks. The problem is that a fixed threshold θ when the size of a network tends to infinity with high probability leads to a complete graph. But real networks are usually sparse.
Therefore, to introduce growing factorization threshold models we use a threshold function θ := θ(n) which depends on the number of nodes n in the network. Then for every value of network size n we have the same parameters except of threshold θ. This means that at every step, when a new node will be added to the graph, some of the existing edges will be removed. In the next section, we will try to find threshold functions which lead to sparse networks.
To preserve readability of the proofs, we consider only the case d = 3 because proofs for higher dimensions can be derived in a similar way. However, we will give not only mean-field approximations but also strict probabilistic proofs, which to the best of our knowledge have not been done for geographical threshold models yet and can be likely applied in the other works too.

Generating sparse networks
The aim of this section is to model sparse growing networks. To do this, we need to find a proper threshold function.
First, we have studied the growth of the real networks. For example, Fig. 1 shows the growth of a citation graph. The data was obtained from the SNAP 1 database. It can be seen that the function y(x) = 4.95x log x − 40x is a good estimation of the growth rate of this network. That is why we decided to focus on the linearithmic or sub-linearithmic growth rate of the model (here and subsequently, by the growth of the model we mean the growth of the number of edges).

Analysis of the expected number of edges
Let M(n) denote the number of edges in the network of size n. To find its expectation, we need the two following lemmas.

Lemma 1
The probability for a node with weight w to be connected to a random node is 1 https://snap.stanford.edu/data/.

Lemma 2 The edge probability in the network is
To improve readability, we moved the proofs of Lemmas 1 and 2 to Appendix. The next theorem shows that our model can have any growth which is less than quadratic.
Theorem 1 Denote as R(n) such function that R(n) = o(n 2 ) and R(n)>0. Then there exists such threshold function θ(n) that the growth of the model is R(n): Proof It easy to check that P e is a continuous function of θ. The intermediate value theorem states that P e (θ) takes any value between P e (θ = 0) = 1/2 and P e (θ = ∞) = 0 at some point within the interval.
Since R(n) = o(n 2 ) and positive, there exists N such that for all n ≥ N, Taking into account Theorem 1, we obtain parameters for the linearithmic and linear growths of the expected number of edges. where A is a constant depending on the Pareto distribution parameters.

Concentration theorem
In this section, we will find the variance of the number of the edges and prove the concentration theorem Proofs of the following lemmas can be found in the Appendix.

Lemma 3
Suppose that x, y and z are random nodes. Let P < be the probability for the node x to be connected to both nodes y and z. Then the variance of the number of edges M is Lemma 4 Suppose that x, y and z are random nodes. Let P < be the probability for the node x to be connected to both nodes y and z. Then Combining these results, we get the following theorem that will be needed to prove the concentration theorem Proof According to Lemmas 3 and 4 in case of θ ≥ w 2 0 , the variance is According to Lemma 2, the expected number of edges is Combining (8) and (6), we obtain Therefore, where C 1 , C 2 , C 3 , A and B are constants depending on the Pareto distribution parameters. Finally, we obtain Theorem 5 Concentration theorem If θ(n) and EM(n) tends to infinity as n → ∞ and where M is the number of edges in the graph.
Proof According to Chebyshev's inequality, we have Let us estimate the right part of the inequality. Using Theorem 4, we get Using the conditions of the theorem, we obtain Combining Theorems 2, 3 and 5, we obtain the following corollary. In this way, we have proved that the number of edges in the graph does not deviate much from its expected value. It means that having the linearithmic or the sub-linearithmic growth of the expected number of edges we also have the same growth for the actual number of edges.

Degree distribution
In this section, we show that our model follows power-law degree distribution with an exponent of 2 and give two proofs. The first is a mean-field approximation. It is usually applied for a fast checking of hypotheses. The second one is a strict probabilistic proof. To the best of our knowledge it has not been considered in the context of the geographic threshold models yet.
To confirm our proofs, we carried out a computer simulation and plotted complementary cumulative distribution of node degree which is shown on Fig. 2. We also used a discrete power-law fitting method, which is described in [2] and implemented in the network analysis package igraph. 2 We obtained α = 2.16, x min = 4 and a quite large p-value of 0.9984 for the Kolmogorov-Smirnov goodness of fit test.

Theorem 6
Let P(k) be the probability of a random node to have a degree k. If n 1 a θ(n) = o(1) , then there exist such constants C 0 and N 0 such that ∀ k(n) : ∀ n > N 0 k(n) < C 0 n we have

Mean-field approximation
This approximation gives power law only for nodes with weights w ≤ θ w 0 . But the expected number of nodes with weights not satisfying this inequality Em is extremely small As it was shown in Lemma 1, the probability of the node � v i = w i � x i with weight w i = w ≤ θ w 0 to have an edge to another random node is Let k i (w) be the degree of the node v i . Then where I stands for the indicator function. As all nodes are independent, we get In the mean-field approximation, we assume that k i (w) is really close to its expectation and we can substitute it by (n − 1)P e (w) in the following expression for the degree distribution P(k) = f (w) dw dk , where f(w) is a density of weights. Thus, Ek i (w) = (n − 1)P e (w).

Fig. 2
Complementary cumulative distribution of node degree n = 3 · 10 5 , � x i ∈ R 3 , w i ∼ Pareto(3, 1), θ = 66.9 Note that we have not used conditions on k(n) and θ(n) yet, they are needed to estimate residual terms in the following rigorous proof.
Proof Degree k i of the node v i is a binomial random variable. Using the probability P e (w) of the node v i with weight w i = w to have an edge to another random node, we can get the probability that k i equals k: To get the total probability, we need to integrate this expression with respect to w Because of P e (w) is a composite function, the integral breaks up into two parts.

Thus,
For estimating I 1 we can use the formula P e (w) = 1 2 w a 0 θ a (a+1) w a from Lemma 1. After making the substitution to integrate with respect to P e (w) and using the incomplete betafunction, we get For I 2 we can derive an upper bound. Note that for w ≥ θ/w 0 we have Therefore, we obtain the following upper estimate (P e (w)) k (1 − P e (w)) n−k−1 aw a 0 w a+1 dw.
We now combine estimates for I 1 , I 2 and the following estimates for the incomplete beta-function:

This gives us
Let us introduce the following notations: Using n θ a (n) = o(1), for k(n) < C 0 n we get . + 1) , If k(n) is a bounded function, then since ε 0 < 1 and ε 1 < 1 we have Since ε x x → 0 for ε < 1 as x → ∞ there exist constants C 0 and N 0 such that for n > N 0 and k(n) < C 0 n we have (ε 1 ) Thus, we obtain Note that regardless of the shape parameter of the Pareto distribution of weights we always generate networks with a degree distribution following a power law with an exponent equals 2. In the next section, we modify our model to change the exponent of the degree destribution and some other properties of the resulting networks.

Model modifications
In this section, we will show how to modify our model to get new properties and how these modifications will affect the degree distribution.

Directed network
Many real networks are directed. To model them and obtain an exponent of the power law that differs from 2, we changed the condition for the existence of an edge. There will be a directed edge (v i , v j ), if and only if As it follows from the next theorem this modification allows us to tune an exponent of the power law.
Theorem 7 Let P out (k) be the probability of an random node to have out-degree k, P in (k) in-degree k. If n max{α,β}/a /θ(n) = o(1), then there exist constants C 0 and N 0 such that ∀k(n) : ∀n > N 0 k(n) < C 0 n we have Proof Here is a proof for the out-degree distribution. The case of the in-degree distribution is similar. First, let us compute P e (w)-the probability of the node � v i = w i � x i with weight w i = w to have an edge to another random node.
Similar to Lemma 1 we get Thus, we obtain Like in Theorem 6, we have The rest of the proof is similar to the corresponding steps of Theorem 6, so we omit details here.
With α = β this model turns into an undirected case with the power-law exponent equals 2 that agrees with Theorem 6.

Functions of dot product
In our model because of the condition w i w j ( � x i , � x j ) ≥ θ ≥ 0 node v i can only be connected to the node v j if an angle between x i and x j is less than π/2. This is a constraint on the possible neighbors of a node that restricts the scope of our model.
We can solve this issue by changing the condition for the existence of an edge: where h : [−1, 1] → R. On Fig. 3 is an example of how it works in R 2 .
increasing function, positive at least in one point from (−1, 1), then there exist constants C 0 and N 0 such that ∀k(n) : ∀n > N 0 k(n) < C 0 n we have

Short scheme of proof
Here is the scheme of proof for the out-degree distribution. The case of the in-degree is similar.
Restrictions on the function h allow us to modify the proof of the directed case. The main difference is a value of the probability P e (w) of a node � v i = w i � x i with the weight w i = w to have an edge to another random node.
To deal with P e (w), we need to compare w 0 with boundaries for each range of θ w α (w ′ ) β 1. If w 0 < θ 1/β w α/β q 1/β , then Last case is w 0 ≥ θ 1/β w α/β r 1/β . But θ(n) grows with n, and for big enough n this inequality will not be satisfied.
It remains only to show that P out (k) = k −2 (1 + o(1)). But now it is easy to see that the influence of every kind of the principal parts of the integral for P e (w) has been already examined in previous theorems for degree distributions. For example, what is proportional to the one we got in Theorem 7. Therefore, we are not giving here additional details.
For example, described class of functions contains functions like e x and x 2m+1 + c, m ∈ N, for a proper constant c.
Of course, not only this small class of functions h(x) has no influence on the degree distribution. For example, it is easy to show that h(x) = x 2m , m ∈ N also has this property. In this way, a proof will be different only in the computation of P e (w).

Conclusion
In our work, we suggest a new model for scale-free networks generation, which is based on the matrix factorization and has a geographical interpretation. We formalize it for fixed size and growing networks. We proof and validate empirically that degree distribution of resulting networks obeys power law with an exponent of 2.
We also consider several extensions of the model. First, we research the case of the directed network and obtain power-law degree distribution with a tunable exponent. Then, we apply different functions to the dot product of latent features vectors, which give us modifications with interesting properties.
Further research could focus on the deep study of latent features vectors distribution. It seems that not only a uniform distribution over the surface of the sphere should be considered because, for example, cities are not uniformly distributed over the surface of Earth. Besides, we want to try other distributions of weights.

Proof of Lemma 3
Let us enumerate pairs of nodes. Each pair of nodes i has an edge indicator I e i . By definition, we have I e 1 , . . . , I e n(n−1)/2 is the sequence of identically distributed random variables, so their expected value is the same and equals to P e . Since EI 2 e i = EI e i = P e , it follows that If edges e i and e j do not have mutual nodes, then I e i and I e j are independent variables. Therefore, E(I e i I e j ) = E(I e i )E(I e j ) = P 2 e . We get EI e(v,w) I e(v,z) is exactly equal to P < .

Proof of Lemma 4
It can be easily seen that .