# Estimating Parameters of a Probabilistic Heterogeneous Block Model via the EM Algorithm

Authors: Marianna Bolla; Ahmed Elbanna
Journal: Journal of Probability and Statistics
 (2015)

---

## Abstract

We introduce a semiparametric block model for graphs, where the within- and between-cluster edge probabilities are not constants within the blocks but are described by logistic type models, reminiscent of the 50-year-old Rasch model and the newly introducedα-β models. Our purpose is to give a partition of the vertices of an observed graph so that the induced subgraphs and bipartite graphs obey these models, where their strongly interlaced parameters give multiscale evaluation of the vertices at the same time. In this way, a profoundly heterogeneous version of the stochastic block model is built via mixtures of the above submodels, while the parameters are estimated with a special EM iteration.


---

## Body

## 1. Introduction

So far many parametric and nonparametric methods have been proposed for so-called community detection in networks. In the nonparametric scenario, hierarchical or spectral methods were applied to maximize the two- or multiway Newman-Girvan modularity [1–4]; more generally, spectral clustering tools (SC), based on Laplacian or modularity spectra, proved to be feasible to find community, anticommunity, or regular structures in networks [5]. In the parametric setup, certain model parameters are estimated, usually via maximizing thelikelihood function of the graph, that is, the joint probability of our observations under the model equations. This so-called ML estimation is a promising method of statistical inference, has solid theoretical foundations [6, 7], and also supports the common-sense goal of accepting parameter values based on which our sample is the most likely.As for the parametric scenario, in the 2010s,α and β models [8, 9] were developed as the unique graph models where thedegree sequence is a sufficient statistic: given the degree sequence, the distribution of the random graph does not depend on the parameters any more (microcanonical distribution over the model graphs). This fact makes it possible to derive the ML estimate of the parameters in a standard way [10]. Indeed, in the context of network data, a lot of information is contained in the degree sequence, though, perhaps, in a more sophisticated way. The vertices may haveclusters (groups or modules) and their memberships may affect their affinity to make ties. We will find groups of the vertices such that the within- and between-cluster edge probabilities admit certain parametric graph models, the parameters of which are highly interlaced. Here the degree sequence is not a sufficient statistic any more, only if it is restricted to the subgraphs. When making inference, we are partly inspired by the stochastic block model, partly by the Rasch model, and by the rectangular analogue of the α-β models.The generalized random graph model, sometimes called stochastic block model (SBM), was first introduced in [11] and discussed later in [12–18]. This model is the generalization of the classical Erdős-Renyi random graph G
n
(
p
), the first random graph of the history introduced in [19] which corresponds to the one-cluster case: between any pair of the n vertices edges come into existence independently, with the same probability p. The graph G
n
(
P
,
P
k
) on n vertices is ageneralized random graph with k
×
k symmetric probability matrix P
=
(
p
u
v
) and proper k-partition P
k
=
(
C
1
,
…
,
C
k
) of the vertices if vertices of C
u and C
v are connected independently, with probability p
u
v
,
1
≤
u
<
v
≤
k; further, any pair of the vertices within C
u is connected with probability p
u
u
(
u
=
1
,
…
,
k
). Therefore, the subgraph of G
n
(
P
,
P
k
) confined to the vertex set C
u is an Erdős-Renyi type random graph, while the bipartite subgraphs connecting vertices of C
u and C
v
(
u
≠
v
) are random bipartite graphs of edge probability p
u
v. Sometimes we refer to P
k asclustering, where C
1
,
…
,
C
k are the clusters. In fact, seminal Szemeredi’s regularity lemma [20] guarantees the existence of such a structure in any huge graph, albeit with an enormously large number of clusters; therefore, it is not applicable for practical purposes.Though in [18] the SBM is called heterogeneous, it is, in fact, a homogeneous one: the probability to make ties is the same within the clusters or between the cluster pairs. Nonetheless, this probability depends on the actual cluster memberships; given the memberships of the vertices, the probability that they are connected is a given constant, for the estimation of which an algorithm (maximizing the likelihood modularity) is proposed in [12] and the the EM algorithm is described in [5]. Here we want to model more complicated within- and between-cluster relations.As for the nonparametric scenario, the role of the degree sequence is the best enhanced in the normalized Laplacian or normalized modularity based SC. In [4] we extended the notion of themodularity matrix to weighted graphs as follows. Let G
=
(
V
,
W
) be anedge-weighted graph on the n-element vertex-set V with the n
×
n symmetric weight matrix W; the entries satisfy w
i
j
=
w
j
i
≥
0, w
i
i
=
0, and they are similarities between the vertex pairs. Themodularity matrix of G is defined as M
=
W
-
d
d
T, where the entries of d are thegeneralized vertex degrees
d
i
=
∑
j
=
1
n
w
i
j
(
i
=
1
,
…
,
n
). Here W is normalized in such a way that ∑
i
=
1
n
∑
j
=
1
n
w
i
j
=
1, an assumption that does not hurt the generality, since the forthcomingnormalized modularity matrix, to be mostly used, is not affected by the scaling of the entries of W:(1)
M
D
=
D
-
1
/
2
M
D
-
1
/
2
,where D
=
diag
⁡
d
1
,
…
,
d
n is the diagonaldegree matrix.In [4], we also introduced the following spectral relaxation technique to approximate the k-partition of the vertices minimizing the within- and between-cluster discrepancies in the spirit of Szemeredi’s regularity lemma. (This discrepancy measures the homogeneity of the clusters; we will not use this notion here; see [21] for more information.) Let the eigenvalues of M
D, enumerated in decreasing absolute values, be 1
>
μ
1
≥
μ
2
≥
⋯
≥
|
μ
n
|
=
0. Assume that μ
k
-
1
>
μ
k, and denote by u
1
,
…
,
u
k
-
1 the corresponding unit-norm, pairwise orthogonal eigenvectors. Let r
1
,
…
,
r
n
∈
R
k
-
1 be the row vectors of the n
×
(
k
-
1
) matrix of column vectors D
-
1
/
2
u
1
,
…
,
D
-
1
/
2
u
k
-
1; they are called k
-
1-dimensional representatives of the vertices. Theweighted
k
-variance of these representatives is defined as (2)
S
~
k
2
=
min
C
1
,
…
,
C
k
⁡
∑
u
=
1
k
∑
j
∈
C
u
d
j
r
j
-
c
u
2
,where c
u
=
1
/
V
o
l
C
u
∑
j
∈
C
u
d
j
r
j is the weighted center of cluster C
u. It is theweighted
k
-means algorithm that gives this minimum, and the point is that the optimum S
~
k is just the minimum distance between the eigensubspace corresponding to μ
0
,
…
,
μ
k
-
1 and the one of the suitably transformed step-vectors over the k-partitions of V. In Chapter 2 of [5] we also discussed that, in view of subspace perturbation theorems, the larger the gap between μ
k
-
1 and μ
k, the smaller S
~
k.Note that the normalized modularity based spectral clustering is the same as the normalized Laplacian based one (see [5]) with the exception that here not only the bottom, but the large absolute value eigenvalues are considered; further, the heterogeneity of the vertex degrees is encoded into the diagonal degree matrix D. With the above technique we embed the vertices into a low dimensional space (spectral relaxation via eigenvectors) and hence perform metric clustering. However, above the spacial location, SC does not assign any parameters to the vertices. Our method to be introduced finds parameters and classifies the vertices at the same time.Here we propose a profoundly heterogeneous block model by carrying on the Rasch model developed more than 50 years ago for evaluating psychological tests [22, 23]. We will call it Logistic Block Model (LBM). Given the number of clusters and a classification of the vertices, we will use the Rasch model for the bipartite subgraphs but the α-β models for the subgraphs themselves and process an iteration (inner cycle) to find the ML estimate of their parameters. Then, based on their contributions to the overall likelihood, we find a new classification of the vertices via taking conditional expectation and using the Bayes rule. Eventually, the two steps are alternated, giving the outer cycle of the iteration.Our algorithm fits into the framework of the EM algorithm [7, 24], in the context of exponential families. The method was originally developed for missing data, and the name comes from the alternatingexpectation (E) andmaximization (M) steps, where in the E-step (assignment phase) we complete the data by substituting for the missing data via taking conditional expectation, while in the M-step (estimation phase) we find the usual ML estimate of the parameters based on the so completed data. The EM algorithm naturally extends to situations, when not the data itself is missing, but it comes from a finite mixture, and the grouping memberships are the missing parameters. This special type of the EM algorithm developed for mixtures is often called collaborative filtering [25, 26] or Gibbs sampling [27], the roots of which method can be traced back to [28].After proving the convergence of the inner cycle to the unique solution of the likelihood equation in each block separately, the convergence of the outer cycle to a local maximum of the likelihood function is easily seen. The advantage of the LBM is that, unlike SC, above clustering the vertices, it also assigns parameters to them, where parameters depend on their cluster memberships. Therefore we call is semiparametric. In the context of social networks, the clusters can be identified with social strata and the parameters with attitudes of people of one group towards people of the other, where attitude is the same for people in the second group but depends on the individual in the first group. The number of clusters is fixed during the iteration, but an initial number of clusters are obtained by SC, via inspecting the normalized modularity spectrum of the graph. We apply the algorithm to randomly generated and real-world data, where the initial clustering was the one obtained by SC. Then we compare the results of SC, SBM, and LBM by the Rand index, and our LBM shows a better agreement with the SC clusters than the SBM. It seems that SC gives a solution close to a local maximum of LBM, which can be regarded as fine tuning of SC. In fact, without good starting clustering, LBM can run into a local maximum (there are many) far from the global one. Therefore, it needs SC starting; however, its advantage is that, in addition, it is able to estimate parameters too.The paper is organized as follows. In Section2 we describe the building blocks of our model. In the context of the α-β models we refer to already proved facts about the existence of the ML estimate, and if it exists, we discuss the algorithm proposed by [9] together with convergence facts; meanwhilewhile, in the context of the β-γ model, we introduce a novel algorithm and prove the convergence of it in the appendix. In Section 3 we use both of the above algorithms for the subgraphs and bipartite subgraphs of our sample graph, and we connect them together in the framework of the EM algorithm. In Section 4, the algorithm is applied to randomly generated and real-world data, while Section 5 is devoted to a brief discussion.


## 2. The Building Blocks of the LBM

Log-linear and logistic type models to describe contingency tables in a combinatorial fashion were proposed, for example, by [11, 29] and widely used in statistics. Together with the Rasch model, they give the foundation of the unweighted graph and bipartite graph models which are the building blocks of our EM iteration.


### 2.1.α
-
β Models for Undirected Random Graphs

With different parameterization, [8, 9] introduced the following random graph model, where the degree sequence is a sufficient statistic. We have an unweighted, undirected random graph on n vertices without loops, such that edges between distinct vertices come into existence independently, but not with the same probability as in the classical Erdős-Renyi model [19]. This random graph can uniquely be characterized by its n
×
n symmetric adjacency matrix A
=
(
A
i
j
) which has zero diagonal and the entries above the main diagonal are independent Bernoulli random variables whose parameters p
i
j
=
P
(
A
i
j
=
1
) obey the following rule. Actually, we formulate this rule for the p
i
j
/
1
-
p
i
j ratios, the so-calledodds: (3)
p
i
j
1
-
p
i
j
=
α
i
α
j
1
≤
i
<
j
≤
n
,where the parameters α
1
,
…
,
α
n are positive reals. This model is called α
model in [9]. With the parameter transformation β
i
=
ln
⁡
α
i
i
=
1
,
…
,
n, it is equivalent to the β
model of [8] which applies to thelogits: (4)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
+
β
j
1
≤
i
<
j
≤
nwith real parameters β
1
,
…
,
β
n.Conversely, the probabilitiesp
i
j and 1
-
p
i
j can be expressed in terms of the parameters, like (5)
p
i
j
=
α
i
α
j
1
+
α
i
α
j
,
1
-
p
i
j
=
1
1
+
α
i
α
jwhose formulas will be intensively used in the subsequent calculations.We are looking for the ML estimate of the parameter vectorα
_
=
(
α
1
,
…
,
α
n
) or β
_
=
(
β
1
,
…
,
β
n
) based on the observed unweighted, undirected graph as a statistical sample. (It may seem that we have a one-element sample here; however, there are n
2 independent random variables, the adjacencies, in the background.)LetD
=
(
D
1
,
…
,
D
n
) denote the degree vector of the above random graph, where D
i
=
∑
j
=
1
n
A
i
j
i
=
1
,
…
,
n. The random vector D, as a function of the sample entries A
i
j’s, is asufficient statistic for the parameter α
_, or, equivalently, for β
_. Roughly speaking, a sufficient statistic itself contains all the information—which can be retrieved from the data—for the parameter. More precisely, a statistic is sufficient when the conditional distribution of the sample, given the statistic, does not depend on the parameter any more. By theNeyman-Fisher factorization theorem [6], a statistic is sufficient if and only if the likelihood function of the sample can be factorized into two parts: one which does not contain the parameter and the other, which includes the parameter, contains the sample entries merely compressed into this sufficient statistic. Consider this factorization of the likelihood function (joint probability of A
i
j’s) in our case. Because of the symmetry of A, this is(6)
L
α
_
A
=
∏
i
=
1
n
-
1
∏
j
=
i
+
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
=
∏
i
=
1
n
∏
j
=
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
1
/
2
=
∏
i
=
1
n
∏
j
=
1
n
p
i
j
1
-
p
i
j
A
i
j
∏
i
=
1
n
∏
j
=
1
n
1
-
p
i
j
1
/
2
=
∏
i
=
1
n
α
i
∑
j
=
1
n
A
i
j
∏
j
=
1
n
α
j
∑
i
=
1
n
A
i
j
∏
i
≠
j
1
-
p
i
j
1
/
2
=
∏
i
≠
j
1
1
+
α
i
α
j
1
/
2
∏
i
=
1
n
α
i
D
i
∏
j
=
1
n
α
j
D
j
1
/
2
=
∏
i
<
j
1
1
+
α
i
α
j
∏
i
=
1
n
α
i
D
i
=
C
α
_
×
∏
i
=
1
n
α
i
D
i
,where we used (5) and the facts that A
i
j
=
A
j
i, p
i
j
=
p
j
i
(
i
<
j
) and A
i
i
=
0, p
i
i
=
0
(
i
=
1
,
…
,
n
). Here the partition function C
α
_
=
∏
i
<
j
1
/
1
+
α
i
α
j only depends on α
_, and the whole likelihood function depends on A
i
j’s merely through D
i’s. Therefore, D is a sufficient statistic. The other factor is constantly 1, indicating that the conditional joint distribution of the entries—given D—is uniform, but we will not make use of this fact. The whole model comes from the so-called log-linear way of model building; see [29]. In [8, 10], the converse statement is also proved: the above α model (reparametrized as β model) is the unique one, where the degree sequence is a sufficient statistic.Let(
a
i
j
) be the matrix of the sample realizations (the adjacency entries of the observed graph), let d
i
=
∑
j
=
1
n
a
i
j be the actual degree of vertex i
(
i
=
1
,
…
,
n
), and let d
=
(
d
1
,
…
,
d
n
) be the observed degree vector. The above factorization also indicates that the joint distribution of the entries belongs to the exponential family, and hence, with natural parameterization [24], the maximum likelihood estimate α
_
^ (or equivalently, β
_
^) is derived from the fact that, with it, the observed degree d
i equals the expected one; that is, E
(
D
i
)
=
∑
i
=
1
n
p
i
j. Therefore, α
_
^ is the solution of the followingmaximum likelihood equation: (7)
d
i
=
∑
j
≠
i
n
α
i
α
j
1
+
α
i
α
j
i
=
1
,
…
,
n
.The ML estimate β
_
^ is easily obtained from α
_
^ via taking the logarithms of its coordinates.Before discussing the solution of the system of (7), let us see what conditions a sequence of nonnegative integers should satisfy so that it could be realized as the degree sequence of a graph. The sequence d
1
,
…
,
d
n of nonnegative integers is calledgraphic if there is an unweighted, undirected graph on n vertices such that its vertex degrees are the numbers d
1
,
…
,
d
n in some order. Without loss of generality, d
i’s can be enumerated in nonincreasing order. The Erdősi-Gallai theorem [30] gives the following necessary and sufficient condition for a sequence to be graphic. The sequence d
1
≥
⋯
≥
d
n
≥
0 of integers is graphic if and only if it satisfies the following two conditions: ∑
i
=
1
n
d
i is even and (8)
∑
i
=
1
k
d
i
≤
k
k
-
1
+
∑
i
=
k
+
1
n
min
⁡
k
,
d
i
,
k
=
1
,
…
,
n
-
1
.Note that for nonnegative (not necessarily integer) real sequences, a continuous analogue of (8) is derived in [8]. For given n, the convex hull of all possible graphic degree sequences is a polytope, to be denoted by D
n. Its extreme points are the so-calledthreshold graphs [31]. It is interesting that for n
=
3 all undirected graphs are threshold, since there are 8 possible graphs on 3 nodes, and there are also 8 vertices of D
3; the n
=
2 case is also not of much interest; therefore we will treat the n
>
3 cases only. The number of vertices of D
n superexponentially grows with n [32]; therefore the problem of characterizing threshold graphs has a high computational complexity. Its facial and cofacial sets are fully described in [10]. Apart from the trivial cases (when there is at least one degree equal to 0 or n
-
1), in [33], the authors give the following equivalent characterization of a threshold graph for n
≥
4: it has no four different vertices, a
,
b
,
c
,
d, such that a
,
b and c
,
d are connected by an edge, but a
,
c and b
,
d are not; that is, it has no two disjoint copies of the complete graph K
2.The authors of [8, 9] prove that D
n is the topological closure of the set of expected degree sequences, and, for given n
>
3, if d
∈
int
⁡
(
D
n
) is an interior point, then maximum likelihood equation (7) has a unique solution. Later, it turned out that the converse is also true: in [10] the authors prove that the ML estimate exists if and only if the observed degree vector is an inner point of D
n. On the contrary, when the observed degree vector is a boundary point of D
n, there is at least one 0 or 1 probability p
i
j which can be obtained only by a parameter vector such that at least one of β
i’s is not finite. In this case, the likelihood function cannot be maximized with a finite parameter set; its supremum is approached with a parameter vector β
_ with at least one coordinate tending to +
∞ or -
∞.The authors in [9] recommend the following algorithm and prove that, provided d
∈
int
⁡
D
n, its iteration converges to the unique solution of system (7). To motivate the iteration, we rewrite (7) as (9)
d
i
=
α
i
∑
j
≠
i
1
1
/
α
j
+
α
i
i
=
1
,
…
,
n
.Then starting with initial parameter values α
1
(
0
)
,
…
,
α
n
(
0
) and using the observed degree sequence d
1
,
…
,
d
n, which is an inner point of D
n, the iteration is as follows:(10)
α
i
t
=
d
i
∑
j
≠
i
1
/
1
/
α
j
t
-
1
+
α
i
t
-
1
i
=
1
,
…
,
nfor t
=
1,2
,
…, until convergence.


### 2.2.β-γ Model for Bipartite Graphs

This bipartite graph model, since there is a one-to-one correspondence between bipartite graphs and 0-1 rectangular arrays, traces back to Haberman [34], Lauritzen [29], and Rasch [22, 23] who applied it for psychological and educational measurements, later market research. According to the Rasch model, the entries of an m
×
n binary table A are independent Bernoulli random variables, where for the parameter p
i
j of the entry A
i
j the following holds: (11)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
-
δ
j
i
=
1
,
…
,
m
;
j
=
1
,
…
,
nwith real parameters β
1
,
…
,
β
m and δ
1
,
…
,
δ
n. As an example, Rasch in [22] investigated binary tables where the rows corresponded to persons and the columns to items of some psychological test, whereas the jth entry of the ith row was 1 if person i answered test item j correctly and 0 otherwise. He also gave a description of the parameters: β
i was the ability of person i, while δ
j was the difficulty of test item j. Therefore, in view of model equation (11), the more intelligent the person is and the less difficult the test is, the larger the success/failure ratio was on a logarithmic scale.Given anm
×
n random binary table A
=
(
A
i
j
), or, equivalently, a bipartite graph, our model is (12)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
+
γ
j
i
=
1
,
…
,
m
,
j
=
1
,
…
,
nwith real parameters β
1
,
…
,
β
m and γ
1
,
…
,
γ
n; further, p
i
j
=
P
(
A
i
j
=
1
).In terms of the transformed parametersb
i
=
e
β
i and g
j
=
e
γ
j, model (12) is equivalent to (13)
p
i
j
1
-
p
i
j
=
b
i
g
j
i
=
1
,
…
,
m
,
j
=
1
,
…
,
n
,where b
1
,
…
,
b
m and g
1
,
…
,
g
n are positive reals.Conversely, the probabilities can be expressed in terms of the parameters:(14)
p
i
j
=
b
i
g
j
1
+
b
i
g
j
,
1
-
p
i
j
=
1
1
+
b
i
g
j
.Observe that if (12) holds with the parameters β
i’s and γ
j’s, then it also holds with the transformed parameters β
i
′
=
β
i
+
c
(
i
=
1
,
…
,
m
) and γ
j
′
=
γ
j
-
c
(
j
=
1
,
…
,
n
) with some c
∈
R. Equivalently, if (13) holds with the positive parameters b
i’s and g
j’s, then it also holds with the transformed parameters (15)
b
i
′
=
b
i
κ
,
g
j
′
=
g
j
κwith some κ
>
0. Therefore, the parameters b
i and g
j are arbitrary to within a multiplicative constant.Here the row-sumsR
i
=
∑
j
=
1
n
A
i
j and the column-sums C
j
=
∑
i
=
1
m
A
i
j are the sufficient statistics for the parameters collected in b
=
b
1
,
…
,
b
m and g
=
g
1
,
…
,
g
n. Indeed, the likelihood function is factorized as(16)
L
b
,
g
A
=
∏
i
=
1
m
∏
j
=
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
=
∏
i
=
1
m
∏
j
=
1
n
p
i
j
1
-
p
i
j
A
i
j
∏
i
=
1
m
∏
j
=
1
n
1
-
p
i
j
=
∏
i
=
1
m
b
i
∑
j
=
1
n
A
i
j
∏
j
=
1
n
g
j
∑
i
=
1
m
A
i
j
∏
i
=
1
m
∏
j
=
1
n
1
-
p
i
j
=
∏
i
=
1
m
∏
j
=
1
n
1
1
+
b
i
g
j
∏
i
=
1
m
b
i
R
i
∏
j
=
1
n
g
j
C
j
.Since the likelihood function depends on A only through its row- and column-sums, by the Neyman-Fisher factorization theorem, R
1
,
…
,
R
m, C
1
,
…
,
C
n are a sufficient statistic for the parameters. The first factor (including the partition function) depends only on the parameters and the row- and column-sums, whereas the seemingly not present factor—which would depend merely on A—is constantly 1, indicating that the conditional joint distribution of the entries, given the row- and column-sums, is uniform (microcanonical) in this model. Note that, in [35], the author characterizes random tables sampled uniformly from the set of 0-1 matrices with fixed margins. Given the margins, the contingency tables coming from the above model are uniformly distributed, and a typical table of this distribution is produced by the β-γ model with parameters estimated via the row- and column-sums as sufficient statistics. In this way, here we obtain another view of the typical table of [35].Based on an observed binary table(
a
i
j
), since we are in exponential family and β
1
,
…
,
β
m, γ
1
,
…
,
γ
n are natural parameters, the likelihood equation is obtained by making the expectation of the sufficient statistic equal to its sample value. Therefore, with the notations r
i
=
∑
j
=
1
n
a
i
j
(
i
=
1
,
…
,
m
) and c
j
=
∑
i
=
1
m
a
i
j
(
j
=
1
,
…
,
n
), the followingsystem of likelihood equations is yielded:(17)
r
i
=
∑
j
=
1
n
b
i
g
j
1
+
b
i
g
j
=
b
i
∑
j
=
1
n
1
1
/
g
j
+
b
i
,
i
=
1
,
…
,
m
;
c
j
=
∑
i
=
1
m
b
i
g
j
1
+
b
i
g
j
=
g
j
∑
i
=
1
m
1
1
/
b
i
+
g
j
,
j
=
1
,
…
,
n
.Note that, for any sample realization of A, (18)
∑
i
=
1
m
r
i
=
∑
j
=
1
n
c
jholds automatically. Therefore, there is dependence between the equations of system (17), indicating that the solution is not unique, in accord with our previous remark about the arbitrary scaling factor κ
>
0 of (15). We will prove that, apart from this scaling, the solution is unique if it exists at all. For our convenience, let (
b
~
,
g
~
) denote the equivalence class of the parameter vector (
b
,
g
), which consists of parameter vectors (
b
′
,
g
′
) satisfying (15) with some κ
>
0. So to avoid this indeterminacy, we may impose conditions on the parameters; for example,(19)
∑
i
=
1
m
β
i
+
∑
j
=
1
n
γ
j
=
0
.Like the graphic sequences, here the following sufficient conditions can be given for the sequencesr
1
≥
⋯
≥
r
m
>
0 and c
1
≥
⋯
≥
c
n
>
0 of integers to be row- and column-sums of an m
×
n matrix of 0-1 entry (see, e.g., [36]):(20)
∑
i
=
1
k
r
i
≤
∑
j
=
1
n
min
⁡
c
j
,
k
,
k
=
1
,
…
,
m
;
∑
j
=
1
k
c
j
≤
∑
i
=
1
m
min
⁡
r
i
,
k
,
k
=
1
,
…
,
n
.Observe that the k
=
1 cases imply r
1
≤
n and c
1
≤
m, whereas the k
=
m and k
=
n cases together imply ∑
i
=
1
m
r
i
=
∑
j
=
1
n
c
j. This statement is the counterpart of the Erdős-Gallai conditions for bipartite graphs, where—due to (18)—the sum of the degrees is automatically even. In fact, the conditions in (20) are redundant: one of the conditions—either the one for the rows or the one for the columns—suffices together with (18) and c
1
≤
m or r
1
≤
n. The so obtained necessary and sufficient conditions definebipartite realizable sequences with the wording of [33]. Already, in 1957, the author of [37] determined arithmetic conditions for the construction of a 0-1 matrix having given row- and column-sums. The construction was given via swaps. More generally, [38] referred to the transportation problem and the Ford-Fulkerson max flow–min cut theorem [39].The convex hull of the bipartite realizable sequencesr
=
(
r
1
,
…
,
r
m
) and c
=
(
c
1
,
…
,
c
n
) forms a polytope in R
m
+
n, actually, because of (18), in its (
m
+
n
-
1
)-dimensional hyperplane. It is calledpolytope of bipartite degree sequences and denoted by P
m
,
n in [33]. It is the special case of the transportation polytope describing margins of contingency tables with nonnegative integer entries. There is an expanding literature on the number of such matrices, for example, [40], and on the number of 0-1 matrices with prescribed row- and column-sums, for example, [41].Analogously to the considerations of theα-β models, and applying the thoughts of the proofs in [8–10], P
m
,
n is the closure of the set of the expected row- and column-sum sequences in the above model. In [33], it is proved that an m
×
n binary table, or equivalently a bipartite graph on the independent sets of m and n vertices, is on the boundary of P
m
,
n if it does not contain two vertex-disjoint edges. In this case, the likelihood function cannot be maximized with a finite parameter set; its supremum is approached with a parameter vector with at least one coordinate, β
i or γ
j, tending to +
∞ or -
∞, or, equivalently, with at least one coordinate, b
i or g
j, tending to +
∞ or 0. Based on the proofs of [10], and stated as Theorem 6.3 in the supplementary material of [10], the maximum likelihood estimate of the parameters of model (13) exists if and only if the observed row- and column-sum sequence r
,
c
∈
r
i
(
P
m
,
n
), the relative interior of P
m
,
n, satisfying (18). In this case for the probabilities, calculated by formula (14) through the estimated positive parameter values b
^
i’s and g
^
j’s (solutions of (17)), 0
<
p
i
j
<
1 holds ∀
i
,
j.Under these conditions, we define an algorithm that converges to the unique (up to the above equivalence) solution of maximum likelihood equations (17).Theorem 1.
If(
r
,
c
)
∈
r
i
(
P
m
,
n
), then the following algorithm gives a unique equivalence class of the parameter vectors as the fixed point of the iteration, which therefore provides the ML estimate of the parameters.
Starting with positive parameter valuesb
i
(
0
)
(
i
=
1
,
…
,
m
) and g
j
(
0
)
(
j
=
1
,
…
,
n
) and using the observed row- and column-sums, the iteration is as follows:(21)
I
b
i
t
=
r
i
∑
j
=
1
n
1
/
1
/
g
j
t
-
1
+
b
i
t
-
1
,
i
=
1
,
…
,
m
I
I
g
j
t
=
c
j
∑
i
=
1
m
1
/
1
/
b
i
t
+
g
j
t
-
1
,
j
=
1
,
…
,
nfor t
=
1,2
,
…, until convergence.


## 2.1.α
-
β Models for Undirected Random Graphs

With different parameterization, [8, 9] introduced the following random graph model, where the degree sequence is a sufficient statistic. We have an unweighted, undirected random graph on n vertices without loops, such that edges between distinct vertices come into existence independently, but not with the same probability as in the classical Erdős-Renyi model [19]. This random graph can uniquely be characterized by its n
×
n symmetric adjacency matrix A
=
(
A
i
j
) which has zero diagonal and the entries above the main diagonal are independent Bernoulli random variables whose parameters p
i
j
=
P
(
A
i
j
=
1
) obey the following rule. Actually, we formulate this rule for the p
i
j
/
1
-
p
i
j ratios, the so-calledodds: (3)
p
i
j
1
-
p
i
j
=
α
i
α
j
1
≤
i
<
j
≤
n
,where the parameters α
1
,
…
,
α
n are positive reals. This model is called α
model in [9]. With the parameter transformation β
i
=
ln
⁡
α
i
i
=
1
,
…
,
n, it is equivalent to the β
model of [8] which applies to thelogits: (4)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
+
β
j
1
≤
i
<
j
≤
nwith real parameters β
1
,
…
,
β
n.Conversely, the probabilitiesp
i
j and 1
-
p
i
j can be expressed in terms of the parameters, like (5)
p
i
j
=
α
i
α
j
1
+
α
i
α
j
,
1
-
p
i
j
=
1
1
+
α
i
α
jwhose formulas will be intensively used in the subsequent calculations.We are looking for the ML estimate of the parameter vectorα
_
=
(
α
1
,
…
,
α
n
) or β
_
=
(
β
1
,
…
,
β
n
) based on the observed unweighted, undirected graph as a statistical sample. (It may seem that we have a one-element sample here; however, there are n
2 independent random variables, the adjacencies, in the background.)LetD
=
(
D
1
,
…
,
D
n
) denote the degree vector of the above random graph, where D
i
=
∑
j
=
1
n
A
i
j
i
=
1
,
…
,
n. The random vector D, as a function of the sample entries A
i
j’s, is asufficient statistic for the parameter α
_, or, equivalently, for β
_. Roughly speaking, a sufficient statistic itself contains all the information—which can be retrieved from the data—for the parameter. More precisely, a statistic is sufficient when the conditional distribution of the sample, given the statistic, does not depend on the parameter any more. By theNeyman-Fisher factorization theorem [6], a statistic is sufficient if and only if the likelihood function of the sample can be factorized into two parts: one which does not contain the parameter and the other, which includes the parameter, contains the sample entries merely compressed into this sufficient statistic. Consider this factorization of the likelihood function (joint probability of A
i
j’s) in our case. Because of the symmetry of A, this is(6)
L
α
_
A
=
∏
i
=
1
n
-
1
∏
j
=
i
+
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
=
∏
i
=
1
n
∏
j
=
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
1
/
2
=
∏
i
=
1
n
∏
j
=
1
n
p
i
j
1
-
p
i
j
A
i
j
∏
i
=
1
n
∏
j
=
1
n
1
-
p
i
j
1
/
2
=
∏
i
=
1
n
α
i
∑
j
=
1
n
A
i
j
∏
j
=
1
n
α
j
∑
i
=
1
n
A
i
j
∏
i
≠
j
1
-
p
i
j
1
/
2
=
∏
i
≠
j
1
1
+
α
i
α
j
1
/
2
∏
i
=
1
n
α
i
D
i
∏
j
=
1
n
α
j
D
j
1
/
2
=
∏
i
<
j
1
1
+
α
i
α
j
∏
i
=
1
n
α
i
D
i
=
C
α
_
×
∏
i
=
1
n
α
i
D
i
,where we used (5) and the facts that A
i
j
=
A
j
i, p
i
j
=
p
j
i
(
i
<
j
) and A
i
i
=
0, p
i
i
=
0
(
i
=
1
,
…
,
n
). Here the partition function C
α
_
=
∏
i
<
j
1
/
1
+
α
i
α
j only depends on α
_, and the whole likelihood function depends on A
i
j’s merely through D
i’s. Therefore, D is a sufficient statistic. The other factor is constantly 1, indicating that the conditional joint distribution of the entries—given D—is uniform, but we will not make use of this fact. The whole model comes from the so-called log-linear way of model building; see [29]. In [8, 10], the converse statement is also proved: the above α model (reparametrized as β model) is the unique one, where the degree sequence is a sufficient statistic.Let(
a
i
j
) be the matrix of the sample realizations (the adjacency entries of the observed graph), let d
i
=
∑
j
=
1
n
a
i
j be the actual degree of vertex i
(
i
=
1
,
…
,
n
), and let d
=
(
d
1
,
…
,
d
n
) be the observed degree vector. The above factorization also indicates that the joint distribution of the entries belongs to the exponential family, and hence, with natural parameterization [24], the maximum likelihood estimate α
_
^ (or equivalently, β
_
^) is derived from the fact that, with it, the observed degree d
i equals the expected one; that is, E
(
D
i
)
=
∑
i
=
1
n
p
i
j. Therefore, α
_
^ is the solution of the followingmaximum likelihood equation: (7)
d
i
=
∑
j
≠
i
n
α
i
α
j
1
+
α
i
α
j
i
=
1
,
…
,
n
.The ML estimate β
_
^ is easily obtained from α
_
^ via taking the logarithms of its coordinates.Before discussing the solution of the system of (7), let us see what conditions a sequence of nonnegative integers should satisfy so that it could be realized as the degree sequence of a graph. The sequence d
1
,
…
,
d
n of nonnegative integers is calledgraphic if there is an unweighted, undirected graph on n vertices such that its vertex degrees are the numbers d
1
,
…
,
d
n in some order. Without loss of generality, d
i’s can be enumerated in nonincreasing order. The Erdősi-Gallai theorem [30] gives the following necessary and sufficient condition for a sequence to be graphic. The sequence d
1
≥
⋯
≥
d
n
≥
0 of integers is graphic if and only if it satisfies the following two conditions: ∑
i
=
1
n
d
i is even and (8)
∑
i
=
1
k
d
i
≤
k
k
-
1
+
∑
i
=
k
+
1
n
min
⁡
k
,
d
i
,
k
=
1
,
…
,
n
-
1
.Note that for nonnegative (not necessarily integer) real sequences, a continuous analogue of (8) is derived in [8]. For given n, the convex hull of all possible graphic degree sequences is a polytope, to be denoted by D
n. Its extreme points are the so-calledthreshold graphs [31]. It is interesting that for n
=
3 all undirected graphs are threshold, since there are 8 possible graphs on 3 nodes, and there are also 8 vertices of D
3; the n
=
2 case is also not of much interest; therefore we will treat the n
>
3 cases only. The number of vertices of D
n superexponentially grows with n [32]; therefore the problem of characterizing threshold graphs has a high computational complexity. Its facial and cofacial sets are fully described in [10]. Apart from the trivial cases (when there is at least one degree equal to 0 or n
-
1), in [33], the authors give the following equivalent characterization of a threshold graph for n
≥
4: it has no four different vertices, a
,
b
,
c
,
d, such that a
,
b and c
,
d are connected by an edge, but a
,
c and b
,
d are not; that is, it has no two disjoint copies of the complete graph K
2.The authors of [8, 9] prove that D
n is the topological closure of the set of expected degree sequences, and, for given n
>
3, if d
∈
int
⁡
(
D
n
) is an interior point, then maximum likelihood equation (7) has a unique solution. Later, it turned out that the converse is also true: in [10] the authors prove that the ML estimate exists if and only if the observed degree vector is an inner point of D
n. On the contrary, when the observed degree vector is a boundary point of D
n, there is at least one 0 or 1 probability p
i
j which can be obtained only by a parameter vector such that at least one of β
i’s is not finite. In this case, the likelihood function cannot be maximized with a finite parameter set; its supremum is approached with a parameter vector β
_ with at least one coordinate tending to +
∞ or -
∞.The authors in [9] recommend the following algorithm and prove that, provided d
∈
int
⁡
D
n, its iteration converges to the unique solution of system (7). To motivate the iteration, we rewrite (7) as (9)
d
i
=
α
i
∑
j
≠
i
1
1
/
α
j
+
α
i
i
=
1
,
…
,
n
.Then starting with initial parameter values α
1
(
0
)
,
…
,
α
n
(
0
) and using the observed degree sequence d
1
,
…
,
d
n, which is an inner point of D
n, the iteration is as follows:(10)
α
i
t
=
d
i
∑
j
≠
i
1
/
1
/
α
j
t
-
1
+
α
i
t
-
1
i
=
1
,
…
,
nfor t
=
1,2
,
…, until convergence.


## 2.2.β-γ Model for Bipartite Graphs

This bipartite graph model, since there is a one-to-one correspondence between bipartite graphs and 0-1 rectangular arrays, traces back to Haberman [34], Lauritzen [29], and Rasch [22, 23] who applied it for psychological and educational measurements, later market research. According to the Rasch model, the entries of an m
×
n binary table A are independent Bernoulli random variables, where for the parameter p
i
j of the entry A
i
j the following holds: (11)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
-
δ
j
i
=
1
,
…
,
m
;
j
=
1
,
…
,
nwith real parameters β
1
,
…
,
β
m and δ
1
,
…
,
δ
n. As an example, Rasch in [22] investigated binary tables where the rows corresponded to persons and the columns to items of some psychological test, whereas the jth entry of the ith row was 1 if person i answered test item j correctly and 0 otherwise. He also gave a description of the parameters: β
i was the ability of person i, while δ
j was the difficulty of test item j. Therefore, in view of model equation (11), the more intelligent the person is and the less difficult the test is, the larger the success/failure ratio was on a logarithmic scale.Given anm
×
n random binary table A
=
(
A
i
j
), or, equivalently, a bipartite graph, our model is (12)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
+
γ
j
i
=
1
,
…
,
m
,
j
=
1
,
…
,
nwith real parameters β
1
,
…
,
β
m and γ
1
,
…
,
γ
n; further, p
i
j
=
P
(
A
i
j
=
1
).In terms of the transformed parametersb
i
=
e
β
i and g
j
=
e
γ
j, model (12) is equivalent to (13)
p
i
j
1
-
p
i
j
=
b
i
g
j
i
=
1
,
…
,
m
,
j
=
1
,
…
,
n
,where b
1
,
…
,
b
m and g
1
,
…
,
g
n are positive reals.Conversely, the probabilities can be expressed in terms of the parameters:(14)
p
i
j
=
b
i
g
j
1
+
b
i
g
j
,
1
-
p
i
j
=
1
1
+
b
i
g
j
.Observe that if (12) holds with the parameters β
i’s and γ
j’s, then it also holds with the transformed parameters β
i
′
=
β
i
+
c
(
i
=
1
,
…
,
m
) and γ
j
′
=
γ
j
-
c
(
j
=
1
,
…
,
n
) with some c
∈
R. Equivalently, if (13) holds with the positive parameters b
i’s and g
j’s, then it also holds with the transformed parameters (15)
b
i
′
=
b
i
κ
,
g
j
′
=
g
j
κwith some κ
>
0. Therefore, the parameters b
i and g
j are arbitrary to within a multiplicative constant.Here the row-sumsR
i
=
∑
j
=
1
n
A
i
j and the column-sums C
j
=
∑
i
=
1
m
A
i
j are the sufficient statistics for the parameters collected in b
=
b
1
,
…
,
b
m and g
=
g
1
,
…
,
g
n. Indeed, the likelihood function is factorized as(16)
L
b
,
g
A
=
∏
i
=
1
m
∏
j
=
1
n
p
i
j
A
i
j
1
-
p
i
j
1
-
A
i
j
=
∏
i
=
1
m
∏
j
=
1
n
p
i
j
1
-
p
i
j
A
i
j
∏
i
=
1
m
∏
j
=
1
n
1
-
p
i
j
=
∏
i
=
1
m
b
i
∑
j
=
1
n
A
i
j
∏
j
=
1
n
g
j
∑
i
=
1
m
A
i
j
∏
i
=
1
m
∏
j
=
1
n
1
-
p
i
j
=
∏
i
=
1
m
∏
j
=
1
n
1
1
+
b
i
g
j
∏
i
=
1
m
b
i
R
i
∏
j
=
1
n
g
j
C
j
.Since the likelihood function depends on A only through its row- and column-sums, by the Neyman-Fisher factorization theorem, R
1
,
…
,
R
m, C
1
,
…
,
C
n are a sufficient statistic for the parameters. The first factor (including the partition function) depends only on the parameters and the row- and column-sums, whereas the seemingly not present factor—which would depend merely on A—is constantly 1, indicating that the conditional joint distribution of the entries, given the row- and column-sums, is uniform (microcanonical) in this model. Note that, in [35], the author characterizes random tables sampled uniformly from the set of 0-1 matrices with fixed margins. Given the margins, the contingency tables coming from the above model are uniformly distributed, and a typical table of this distribution is produced by the β-γ model with parameters estimated via the row- and column-sums as sufficient statistics. In this way, here we obtain another view of the typical table of [35].Based on an observed binary table(
a
i
j
), since we are in exponential family and β
1
,
…
,
β
m, γ
1
,
…
,
γ
n are natural parameters, the likelihood equation is obtained by making the expectation of the sufficient statistic equal to its sample value. Therefore, with the notations r
i
=
∑
j
=
1
n
a
i
j
(
i
=
1
,
…
,
m
) and c
j
=
∑
i
=
1
m
a
i
j
(
j
=
1
,
…
,
n
), the followingsystem of likelihood equations is yielded:(17)
r
i
=
∑
j
=
1
n
b
i
g
j
1
+
b
i
g
j
=
b
i
∑
j
=
1
n
1
1
/
g
j
+
b
i
,
i
=
1
,
…
,
m
;
c
j
=
∑
i
=
1
m
b
i
g
j
1
+
b
i
g
j
=
g
j
∑
i
=
1
m
1
1
/
b
i
+
g
j
,
j
=
1
,
…
,
n
.Note that, for any sample realization of A, (18)
∑
i
=
1
m
r
i
=
∑
j
=
1
n
c
jholds automatically. Therefore, there is dependence between the equations of system (17), indicating that the solution is not unique, in accord with our previous remark about the arbitrary scaling factor κ
>
0 of (15). We will prove that, apart from this scaling, the solution is unique if it exists at all. For our convenience, let (
b
~
,
g
~
) denote the equivalence class of the parameter vector (
b
,
g
), which consists of parameter vectors (
b
′
,
g
′
) satisfying (15) with some κ
>
0. So to avoid this indeterminacy, we may impose conditions on the parameters; for example,(19)
∑
i
=
1
m
β
i
+
∑
j
=
1
n
γ
j
=
0
.Like the graphic sequences, here the following sufficient conditions can be given for the sequencesr
1
≥
⋯
≥
r
m
>
0 and c
1
≥
⋯
≥
c
n
>
0 of integers to be row- and column-sums of an m
×
n matrix of 0-1 entry (see, e.g., [36]):(20)
∑
i
=
1
k
r
i
≤
∑
j
=
1
n
min
⁡
c
j
,
k
,
k
=
1
,
…
,
m
;
∑
j
=
1
k
c
j
≤
∑
i
=
1
m
min
⁡
r
i
,
k
,
k
=
1
,
…
,
n
.Observe that the k
=
1 cases imply r
1
≤
n and c
1
≤
m, whereas the k
=
m and k
=
n cases together imply ∑
i
=
1
m
r
i
=
∑
j
=
1
n
c
j. This statement is the counterpart of the Erdős-Gallai conditions for bipartite graphs, where—due to (18)—the sum of the degrees is automatically even. In fact, the conditions in (20) are redundant: one of the conditions—either the one for the rows or the one for the columns—suffices together with (18) and c
1
≤
m or r
1
≤
n. The so obtained necessary and sufficient conditions definebipartite realizable sequences with the wording of [33]. Already, in 1957, the author of [37] determined arithmetic conditions for the construction of a 0-1 matrix having given row- and column-sums. The construction was given via swaps. More generally, [38] referred to the transportation problem and the Ford-Fulkerson max flow–min cut theorem [39].The convex hull of the bipartite realizable sequencesr
=
(
r
1
,
…
,
r
m
) and c
=
(
c
1
,
…
,
c
n
) forms a polytope in R
m
+
n, actually, because of (18), in its (
m
+
n
-
1
)-dimensional hyperplane. It is calledpolytope of bipartite degree sequences and denoted by P
m
,
n in [33]. It is the special case of the transportation polytope describing margins of contingency tables with nonnegative integer entries. There is an expanding literature on the number of such matrices, for example, [40], and on the number of 0-1 matrices with prescribed row- and column-sums, for example, [41].Analogously to the considerations of theα-β models, and applying the thoughts of the proofs in [8–10], P
m
,
n is the closure of the set of the expected row- and column-sum sequences in the above model. In [33], it is proved that an m
×
n binary table, or equivalently a bipartite graph on the independent sets of m and n vertices, is on the boundary of P
m
,
n if it does not contain two vertex-disjoint edges. In this case, the likelihood function cannot be maximized with a finite parameter set; its supremum is approached with a parameter vector with at least one coordinate, β
i or γ
j, tending to +
∞ or -
∞, or, equivalently, with at least one coordinate, b
i or g
j, tending to +
∞ or 0. Based on the proofs of [10], and stated as Theorem 6.3 in the supplementary material of [10], the maximum likelihood estimate of the parameters of model (13) exists if and only if the observed row- and column-sum sequence r
,
c
∈
r
i
(
P
m
,
n
), the relative interior of P
m
,
n, satisfying (18). In this case for the probabilities, calculated by formula (14) through the estimated positive parameter values b
^
i’s and g
^
j’s (solutions of (17)), 0
<
p
i
j
<
1 holds ∀
i
,
j.Under these conditions, we define an algorithm that converges to the unique (up to the above equivalence) solution of maximum likelihood equations (17).Theorem 1.
If(
r
,
c
)
∈
r
i
(
P
m
,
n
), then the following algorithm gives a unique equivalence class of the parameter vectors as the fixed point of the iteration, which therefore provides the ML estimate of the parameters.
Starting with positive parameter valuesb
i
(
0
)
(
i
=
1
,
…
,
m
) and g
j
(
0
)
(
j
=
1
,
…
,
n
) and using the observed row- and column-sums, the iteration is as follows:(21)
I
b
i
t
=
r
i
∑
j
=
1
n
1
/
1
/
g
j
t
-
1
+
b
i
t
-
1
,
i
=
1
,
…
,
m
I
I
g
j
t
=
c
j
∑
i
=
1
m
1
/
1
/
b
i
t
+
g
j
t
-
1
,
j
=
1
,
…
,
nfor t
=
1,2
,
…, until convergence.


## 3. Parameter Estimation in the LBM

In the several clusters’ case, we are putting the bricks together. The above discussedα-β and β-γ models will be the building blocks of the LBM to be introduced. Here the degree sequences are not any more sufficient for the whole graph, only for the building blocks of the subgraphs.Given1
≤
k
≤
n, we are looking for k-partition, in other words, clusters C
1
,
…
,
C
k of the vertices, such that(i)
different vertices are independently assigned to a clusterC
u with probability π
u
(
u
=
1
,
…
,
k
), where ∑
u
=
1
k
π
u
=
1;
(ii)
given the cluster memberships verticesi
∈
C
u and j
∈
C
v are connected independently, with probability p
i
j such that (22)
ln
⁡
p
i
j
1
-
p
i
j
=
β
i
v
+
β
j
u
,
for any1
≤
u, v
≤
k pair. Equivalently, (23)
p
i
j
1
-
p
i
j
=
b
i
c
j
b
j
c
i
,
wherec
i is the cluster membership of vertex i and b
i
v
=
e
β
i
v.The parameters are collected in the vectorπ
_
=
(
π
1
,
…
,
π
k
) and the n
×
k matrix B of b
i
u’s (
i
∈
C
u
,
u
=
1
,
…
,
k
). The likelihood function is the following mixture: (24)
∑
1
≤
u
,
v
≤
k
π
u
π
v
∏
i
∈
C
u
,
j
∈
C
v
p
i
j
a
i
j
1
-
p
i
j
1
-
a
i
j
.Here A
=
(
a
i
j
) is the incomplete data specification as the cluster memberships are missing. Therefore, it is straightforward to use the EM algorithm, proposed by [24], also discussed in [7, 42], for parameter estimation from incomplete data. This special application for mixtures is sometimes calledcollaborative filtering (see [25, 26]) which is rather applicable to fuzzy clustering.First we complete our data matrixA with latent membership vectors m
1
,
…
,
m
n of the vertices that are k-dimensional i.i.d. Multy
(
1
,
π
_
) (multinomially distributed) random vectors. More precisely, m
i
=
(
m
i
1
,
…
,
m
i
k
), where m
i
u
=
1 if i
∈
C
u and zero otherwise. Thus, the sum of the coordinates of any m
i is 1, and P
(
m
i
u
=
1
)
=
π
u.Note that if the cluster memberships were known, then the complete likelihood would be(25)
∏
u
=
1
k
∏
i
=
1
n
∏
v
=
1
k
∏
j
=
1
n
p
i
j
m
j
v
a
i
j
·
1
-
p
i
j
m
j
v
1
-
a
i
j
m
i
uthat is valid only in case of known cluster memberships.Starting with initial parameter valuesπ
_
(
0
) and B
(
0
) and membership vectors m
1
(
0
)
,
…
,
m
n
(
0
), the tth step of the iteration is the following (t
=
1,2
,
…):(i)
E-step: we calculate the conditional expectation of eachm
i conditioned on the model parameters and on the other cluster assignments obtained in step t
-
1 and collectively denoted by M
(
t
-
1
).
The responsibility of vertexi for cluster u in the tth step is defined as the conditional expectation π
i
u
(
t
)
=
E
m
i
u
∣
M
(
t
-
1
), and, by the Bayes theorem, it is (26)
π
i
u
t
=
P
M
t
-
1
∣
m
i
u
=
1
·
π
u
t
-
1
∑
v
=
1
k
P
M
t
-
1
∣
m
i
v
=
1
·
π
v
t
-
1
(
u
=
1
,
…
,
k
;
i
=
1
,
…
,
n
). For each i, π
i
u
(
t
) is proportional to the numerator; therefore the conditional probabilities P
(
M
(
t
-
1
)
|
m
i
u
=
1
) should be calculated for u
=
1
,
…
,
k. But this is just the part of likelihood (25) affecting vertex i under the condition m
i
u
=
1. Therefore, (27)
P
M
t
-
1
∣
m
i
u
=
1
=
∏
v
=
1
k
∏
j
∈
C
v
,
j
~
i
b
i
v
t
-
1
b
j
u
t
-
1
1
+
b
i
v
t
-
1
b
j
u
t
-
1
∏
j
∈
C
v
,
j
≁
i
1
1
+
b
i
v
t
-
1
b
j
u
t
-
1
=
∏
v
=
1
k
b
i
v
t
-
1
b
j
u
t
-
1
1
+
b
i
v
t
-
1
b
j
u
t
-
1
e
v
i
1
1
+
b
i
v
t
-
1
b
j
u
t
-
1
C
v
·
C
v
-
1
/
2
-
e
v
i
,
wheree
v
i is the number of edges within C
v that are connected to i.
(ii)
M-step: we updateπ
_
(
t
) and m
(
t
): π
u
(
t
)
≔
1
/
n
∑
i
=
1
n
π
i
u
(
t
) and m
i
u
(
t
)
=
1 if π
i
u
(
t
)
=
m
a
x
v
π
i
v
(
t
) and 0 otherwise (in case of ambiguity, we select the smallest index for the cluster membership of vertex i). This is an ML estimation (discrete one, in the latter case, for the cluster membership). In this way, new clustering of the vertices is obtained.
Then we estimate the parameters in the actual clustering of the vertices. In the within-cluster scenario, we use the parameter estimation of model (3), obtaining estimates of b
i
u’s (i
∈
C
u) in each cluster separately (
u
=
1
,
…
,
k
); as for cluster u, b
i
u corresponds to α
i and the number of vertices is C
u. In the between-cluster scenario, we use bipartite graph model (13) in the following way. For u
<
v, edges connecting vertices of C
u and C
v form a bipartite graph, based on which the parameters b
i
v
i
∈
C
u and b
j
u
j
∈
C
v are estimated with the above algorithm; here b
i
v’s correspond to b
i’s, b
j
u’s correspond to g
j’s, and the number of rows and columns of the rectangular array corresponding to this bipartite subgraph of A is C
u and C
v, respectively. With the estimated parameters, collected in the n
×
k matrix B
(
t
), we go back to the E-step and so forth.As in the M-step we increase the likelihood in all parts, and in the E-step we relocate the vertices into the cluster where their likelihoods are maximized, the nonnegative likelihood function is increased in each iteration. Since the likelihood function is bounded from above (unless in some inner cycle we start from the boundary of a polytope of Section2), it must converge to a local maximum.Note that here the parameterβ
i
v with c
i
=
u embodies the affinity of vertex i of cluster C
u towards vertices of cluster C
v; and likewise, β
j
u with c
j
=
v embodies the affinity of vertex j of cluster C
v towards vertices of cluster C
u. By the model, these affinities are added together on the level of the logits. This so-called k-β model, introduced in [43], is applicable to social networks, where attitudes of individuals in the same social group (say, u) are the same toward members of another social group (say, v), though this attitude also depends on the individual in group u. The model may also be applied to biological networks, where the clusters correspond, for example, to different functioning synopses or other units of the brain; see [44].After normalizing theβ
i
v
(
i
∈
C
u
) and β
j
u
(
j
∈
C
v
) to meet the requirement of (19), for any u
≠
v pair, the sum of the parameters will be zero, and their sign and magnitude indicate the affinity of nodes of C
u to make ties with the nodes of C
v and vice versa: (28)
∑
i
∈
C
u
β
i
v
+
∑
j
∈
C
v
β
j
u
=
0
.This becomes important when we want to compare the parameters corresponding to different cluster pairs. For selecting the initial number of clusters, we can use considerations of [45], while for the initial clustering, we can use spectral clustering tools of [5].


## 4. Applications

Now we illustrate the performance of our algorithm via randomly generated and real-world data. Note that while processing the iteration, we sometimes run into threshold subgraphs or bipartite subgraphs on the boundary of the polytope of bipartite degree sequences. Even in this case, our iteration converged for most coordinates of the parameter vectors, while someb
i
v coordinates tended to +
∞ or 0 (numerically, when stopping the iteration, they took on a very “large” or “small” value). This means that the affinity of node i towards nodes of the cluster j is infinitely “large” or “small”; that is, this node is liable to always or never make ties with nodes of cluster j.First we generated a random graph onn
=
580 vertices and with k
=
3 underlying vertex-clusters C
1, C
2, C
3 in the following way. Let C
1
≔
190, C
2
≔
193, and C
3
≔
197. The parameters β
i
1
(
i
∈
C
1
), β
i
2
(
i
∈
C
1
), and β
i
3
(
i
∈
C
1
) were chosen independently at uniform from the intervals [
0,1
], [
-
1,1
], and [
-
1,0.5
], respectively. The parameters β
i
1
(
i
∈
C
2
), β
i
2
(
i
∈
C
2
), and β
i
3
(
i
∈
C
2
) were chosen independently at uniform from the intervals [
-
0.75,0.5
], [
-
1,0
], and [
-
0.5,1
], respectively. The parameters β
i
1
(
i
∈
C
3
), β
i
2
(
i
∈
C
3
), and β
i
3
(
i
∈
C
3
) were chosen independently at uniform from the intervals [
-
0.25,0.75
], [
-
0.25,0.25
], and [
-
0.5,0.5
], respectively.Starting with 3 clusters, obtained by spectral clustering tools, and initial parameter values collected inB
(
0
) of all 1 entries, after some outer steps, the iteration converged to B
^
=
(
b
^
i
v
). With β
^
i
v
=
ln
⁡
b
^
i
v, we plotted the β
i
v, β
^
i
v pairs for i
∈
C
u, u
,
v
=
1,2
,
3. Figure 1 shows a good fit with M
S
E
=
1.14634 of the estimated parameters to the original ones. Indeed, by the general theory of the ML estimation [6], for “large” n, the ML estimate should approach the true parameter, based on which the model was generated. So the good fit means that our algorithm finds estimates close to the true parameters.Figure 1
Data were generated based on parametersβ
i
v’s chosen uniformly in different intervals, k
=
3, C
1
=
190, C
2
=
193, and C
3
=
197. The estimated versus the original parameters β
i
v’s are shown for i
∈
C
u (u
,
v
=
1
,
…
,
k), where β
i
1
~
U
[
0,1
] (i
∈
C
1), β
i
1
~
U
-
0.75,0.5 (i
∈
C
2), β
i
1
~
U
-
0.25,0.75 (i
∈
C
3), β
i
2
~
U
[
-
1,1
] (i
∈
C
1), β
i
2
~
U
[
-
1,0
] (i
∈
C
2), β
i
2
~
U
[
-
0.25,0.25
] (i
∈
C
3), β
i
3
~
U
[
-
1,0.5
] (i
∈
C
1), β
i
3
~
U
[
-
0.5,1
] (i
∈
C
2), and β
i
3
~
U
[
-
0.5,0.5
] (i
∈
C
3), respectively. MSE = 1.14634.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Then we generated the same size of a random graph, where the initial parameters followed Gaussian distribution, with different parameters for the within- and between-cluster relations. Based on the parameters we calculated the edge probabilities, and we generated a random graph with them. Eventually, we estimated the parameters with our algorithm; see Figure2. The Gaussian data are about in the same ranges as the uniform ones; however, they are better concentrated to their means. It can be the cause of a bit smaller M
S
E
=
1.12556.Figure 2
Data were generated based on parametersβ
i
v’s following Gaussian distribution with different parameters for the within- and between-cluster relations, k
=
3, C
1
=
190, |
C
2
|
=
193, and |
C
3
|
=
197. The estimated versus the original parameters β
i
v’s are shown for i
∈
C
u (u
,
v
=
1
,
…
,
k), where β
i
1
~
N
(
0.5,0.25
) (i
∈
C
1), β
i
1
~
N
(
-
0.125,0.312
) (i
∈
C
2), β
i
1
~
N
(
0.25,0.25
) (i
∈
C
3), β
i
2
~
N
(
0,0.5
) (i
∈
C
1), β
i
2
~
N
(
-
0.5,0.25
) (i
∈
C
2), β
i
2
~
N
(
0,0.125
) (i
∈
C
3), β
i
3
~
N
(
-
0.25,0.375
) (i
∈
C
1), β
i
3
~
N
(
0.25,0.375
) (i
∈
C
2), and β
i
3
~
N
(
0,0.25
) (i
∈
C
3), respectively. MSE = 1.12556.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure3 shows the resulting clusters obtained by applying the LBM algorithm to the B&K fraternity data [46] with n
=
58 vertices; see also http://vlado.fmf.uni-lj.si/pub/networks/data/ucinet/ucidata.htm#bkfrat. The data, collected by Bernard and Killworth, are behavioral frequency counts, based on communication frequencies between students of a college fraternity in Morgantown, West Virginia. We used the binarized version of the symmetric edge-weight matrix. When the data were collected, the 58 occupants had been living together for at least three months, but senior students had been living there for up to three years. We used our normalized modularity based spectral clustering algorithm [4] to find the starting clusters. In the normalized modularity spectrum we found a gap after the third eigenvalue (in decreasing order of their absolute values); therefore we applied the algorithm with k
=
4 clusters. The four groups are likely to consist of persons living together for about the same time period.Figure 3
The 4 clusters found by the LBM algorithm in the B&K fraternity data, with 10, 9, 20, and 19 students in the clusters, respectively. RAND index = 1 between the SC and LBM.While processing the iteration, occasionally we bumped into the situation when the degree sequence lied on the boundary of the convex polytopes defined in Sections2.1 and 2.2. Unfortunately, this can occur when our graph is not dense enough. In these situations, the iteration did not converge for some coordinates β
i
v, but they seemed to tend to +
∞ or -
∞. Equivalently, the corresponding b
i
v for some i
∈
C
u and v tended to +
∞ or 0, yielding the situation that member i
∈
C
u had +
∞ or 0 affinity towards members of C
v. Another way, recommended in [8], is to add a small amount to each degree to avoid this situation. However, we did not want to manipulate the original graph, which was too sparse to produce degree sequences in the interior of one or more polytopes.We also analyzed the network based on the friendships between the users of theLast.fm music recommendation system [47]. Last.fm is an online service in music based social networking. Each user may have friends inside the Last.fm social network, and so they form a timestamped undirected graph. In 2012, there were 71,000 users and 285,241 edges between them. Actually, we only used the 15-core of this graph for clustering. Starting with SC, the LBM iteration refined the three underlying clusters; see Figure 4.Figure 4
The 3 clusters found by the LBM algorithm in the network of theLast.fm users with 1012, 97, and 53 users in the clusters, respectively. RAND index = 0.99205 between the SC and LBM.We clustered the vertices of the above networks with the EM algorithm and estimated the parameters in both the LBM and SBM (for the iteration of the letter one, see Chapter 5 of [5]). For the initial clustering we used SC with number of clusters k such that we found a remarkable gap in the normalized modularity spectrum between |
μ
k
-
1
| and |
μ
k
|. During the iteration some clusters may become empty which reduces k; it was not the case in our iterations. It is also possible to start with different values of k; here we started with the smallest possible k which indicated a gap. In case of the B&K fraternity data, the leading eigenvalues in decreasing absolute values (apart from the trivial 1) were μ
1
=
0.235826, μ
2
=
-
0.228652, μ
3
=
0.223039, μ
4
=
-
0.198867, and μ
5
=
0.194783, indicating a gap between |
μ
3
| and |
μ
4
|, so we selected k
=
4. In case of the Last.fm data, the leading eigenvalues in decreasing absolute values (apart from the trivial 1) were μ
1
=
0.97061, μ
2
=
0.942929, μ
3
=
0.892111, and μ
4
=
0.862594, indicating a gap between μ
2 and μ
3, so we selected k
=
3.After some outer iterations both the LBM and SBM converged to a local maximum. We compared the clustering obtained by SC versus LBM and SC versus SBM via the Rand index introduced in [48]. This index is between 0 and 1, and the larger it is, the better the agreement between the two clusterings is. We found a good agreement between the SC clusters and those of the LBM: R
A
N
D
=
1 in Figure 3 and R
A
N
D
=
0.99205 in Figure 4, whereas, between the SC clustering and SBM clustering, we obtained R
A
N
D
=
0.61525 for the B&K fraternity data and R
A
N
D
=
0.96912 for the Last.fm data. This shows that LBM is better fine tuning of the spectral clustering than SBM, at least, in these examples, where the diversity of the degrees is present.


## 5. Discussion

Our model is a profoundly heterogeneous kind of a block model, where the subgraphs and bipartite subgraphs obey parametric graph models, within which the connections are mainly determined by the degrees. The EM type algorithm introduced here finds the blocks and estimates the parameters at the same time.When investigating controllability of large networks, the authors of [49] observe and prove that a system’s controllability is to a great extent encoded by the underlying network’s degree distribution. In our model, this is true only for the building blocks. Possibly, the blocks could be controlled separately, based on the degree sequences of the subgraphs.Our model is applicable to large inhomogeneous networks, and above finding clusters of the vertices, it also assigns multiscale parameters to them. In social networks, these parameters can be associated with attitudes of persons of one group towards those in the same or another group. The attitudes are, in fact, affinities to make ties.We prove the convergence of the inner cycle of the algorithm to the unique solution of the maximum likelihood equation within the subgraphs and bipartite subgraphs. Then, by the iteration of the EM algorithm, the convergence of the outer cycle to a local maximum of the overall likelihood is straightforward. As there can be several local maxima, good starting is important. Our final clusters show a good agreement with the spectral clusters; therefore, the algorithm can be considered as a refinement of the spectral clustering and gives estimates of the parameters which provide a local maximum of the overall LBM likelihood with clusters near to the spectral ones.



---

*Source: 657965-2015-12-24.xml*