COMMUNITY DETECTION WITH NODAL INFORMATION
B Y H AOLEI W ENG AND YANG F ENG∗

arXiv:1610.09735v2 [stat.ME] 10 Dec 2016

Columbia University
Community detection is one of the fundamental problems in the study
of network data. Most existing community detection approaches only consider edge information as inputs, and the output could be suboptimal when
nodal information is available. In such cases, it is desirable to leverage nodal
information for the improvement of community detection accuracy. Towards
this goal, we propose a flexible network model incorporating nodal information, and develop likelihood-based inference methods. For the proposed
methods, we establish favorable asymptotic properties as well as efficient algorithms for computation. Numerical experiments show the effectiveness of
our methods in utilizing nodal information across a variety of simulated and
real network data sets.

1. Introduction. Networked systems are ubiquitous in modern society. Examples include worldwide web, gene regulatory networks, and social networks.
Network analysis has attracted a lot of research attention from social science,
physics, computer science and mathematical science. There have been some interesting findings regarding the network structures, such as small world phenomena and power-law degree distributions (Newman, 2003). One of the fundamental
problems in network analysis is detecting and characterizing community structure
in networks. Communities can be intuitively understood as groups of nodes which
are densely connected within groups while sparsely connected between groups1 .
Identifying network communities not only helps better understand structural fea∗

Partially supported by NSF CAREER Grant DMS-1554804.
AMS 2000 subject classifications: Primary 62F99; secondary 62P25.

Keywords and phrases: Networks, Community detection, Stochastic block model, Multi-logistic
regression, Maximum likelihood, Profile likelihood, Variational inference, Consistency, Semidefinite
programming.
1
More rarely, one can encounter communities of the opposite meaning in disassortative mixing
networks.

1

2

H. WENG AND Y. FENG

tures of the network, but also offers practical benefits. For example, communities
in social networks tend to share similar interest, which could provide useful information to build recommendation systems.
Existing community detection methods can be roughly divided into algorithmic
and model-based ones (Zhao et al., 2012). Algorithmic methods typically define
an objective function such as modularity (Newman, 2006), which measures the
goodness of a network partition, and design algorithms to search for the solution of
the corresponding optimization problem. See Fortunato (2010) for a thorough discussion of various algorithms. Unlike algorithmic approaches, model-based methods first construct statistical models that are assumed to generate the networks
under study, and then develop statistical inference tools to learn the latent communities. Some popular models include stochastic block model (Holland et al., 1983),
degree-corrected stochastic block model (Dasgupta et al., 2004; Karrer and Newman, 2011) and mixed membership stochastic block model (Airoldi et al., 2009).
In recent years, there have been increasingly active researches towards understanding the theoretical performances of community detection methods under different types of models. Regarding the stochastic block model, consistency results
have been proved for likelihood based approaches, including maximum likelihood
(Celisse et al., 2012; Choi et al., 2012), profile likelihood (Bickel and Chen, 2009),
pseudo likelihood (Amini et al., 2013) and variational inference (Celisse et al.,
2012; Bickel et al., 2013), among others. Some of the existing results are generalized to degree-corrected block models (Zhao et al., 2012). Another line of theoretical works focuses on methods of moments. See Rohe et al. (2011); Lei and
Rinaldo (2014); Jin (2015); Qin and Rohe (2013); Joseph and Yu (2013) on theoretical analysis of spectral clustering for detecting communities in block models.
Spectral clustering (Zhang et al., 2014) and tensor spectral method (Anandkumar
et al., 2014) have also been used to detect overlapping communities under mixed
membership models. In addition, carefully constructed convex programming has
been shown to enjoy provable guarantees for community detection (Chen et al.,
2012; Amini and Levina, 2014; Cai and Li, 2015; Guédon and Vershynin, 2016;
Chen et al., 2015). See also the interesting theoretical works of community detec-

COMMUNITY DETECTION WITH NODAL INFORMATION

3

tion under minimax framework (Zhang and Zhou, 2015; Gao et al., 2015). Finally,
there exists a different research theme focusing on detectability instead of consistency (Decelle et al., 2011; Krzakala et al., 2013; Saade et al., 2014; Abbe and
Sandon, 2015).
All the aforementioned methods are based on only the observations of the edge
connections in the networks. In the real world, however, networks often appear
with additional nodal information. For example, social networks such as Facebook
and Twitter contain users’ personal profile information. A citation network has the
authors’ names, keywords, and abstracts of papers. Since nodes in the same communities tend to share similar features, we can expect that nodal attributes are
in turn indicative of community structures. Combining both sources of edge and
nodal information opens the possibility for more accurate community discovery.
Many efficient heuristic algorithms are proposed in recent years to accomplish this
goal (Akoglu et al., 2012; Ruan et al., 2013; Chang and Blei, 2010; Nallapati and
Cohen, 2008; Yang et al., 2013). However, not much theory has been established
to understand the statistical properties. See Binkiewicz et al. (2014); Zhang et al.
(2013) for some theoretical developments. In this paper, we aim to give a thorough
study of the community detection with nodal information problem. Our work first
introduces a flexible modeling framework tuned for community detection when
edge and nodal information coexist. Under a specific model, we then study three
likelihood methods and derive their asymptotic properties. Regarding the computation of the estimators, we resort to a convex relaxation (semidefinite programming)
approach to obtain a preliminary community estimate serving as a good initialization fed into “coordinate” ascent type iterative algorithms, to help locate the
global optima. Various numerical experiments demonstrate that our methods can
accurately discover community structures by making efficient use of nodal information.
The rest of the paper is organized as follows. Section 2 introduces our network
model with nodal information. We then propose likelihood based methods and derive the corresponding asymptotic properties in Section 3. Section 4 is devoted to
the design and analysis of practical algorithms. Simulation examples and real data

4

H. WENG AND Y. FENG

analysis are presented in Section 5. We conclude the paper with a discussion in
Section 6. All the technical proofs are collected in the Appendix.
2. Network Modeling with Nodal Information. A network is usually represented by a graph G(V, E), where V = {1, 2, . . . , n} is the set of nodes and E
is the set of edges. Throughout the paper, we will focus on the networks in which
the corresponding graphs are undirected and contain no self-edges. The observed
edge information can be recorded in the adjacency matrix A ∈ {0, 1}n×n , where
Aij = Aji = 1 if and only if (i, j) ∈ E. Suppose the network can be divided
into K non-overlapping communities. Let c = (c1 , . . . , cn ) be the community assignment vector, with ci denoting the community membership of node i and taking
values in {1, 2, . . . , K}. Additionally, the available nodal information is formulated in a covariate matrix X = (x1 , . . . , xn )T ∈ Rn×p , where xi ∈ Rp is the i-th
node’s covariate vector. The goal is to estimate c from the observations A and X.
2.1. Conditional Independence. We treat A, X and c as random and posit a
statistical model for them. Before introducing the model, we would like to elucidate the main motivation. For the purpose of community detection, we follow the
standard two-step inference procedure:
(1) Derive parameter estimator θ̂ based on P (A, X; θ),
(2) Perform posterior inference according to P (c | A, X; θ̂).
Under this framework, we now make a conditional independence assumption: A ⊥
X | c. Admittedly, the assumption imposes a strong constraint that given the community membership, what nodes are like (described by covariates X) does not
affect how they are connected (encoded in A). On the other hand, this assumption is consistent with our belief that knowing nodal information can help identify
community structure c. More importantly, this assumed conditional independence
turns out to simplify the above two steps to a great extent. First, for the parameter
estimation step, the conditional independence assumption implies that
(1)
P (A, X; θ) =

X
c

P (A | c)P (X | c)P (c) = P (X; θ1 )

X
c

P (A | c; θ2 )P (c | X; θ3 ),

COMMUNITY DETECTION WITH NODAL INFORMATION

5

where θ = (θ1 , θ2 , θ3 ) indexes a family of generative models (not restricted to
parametric forms). Regarding the second step, conditional independence leads to
(2)
P (A | c)P (c | X)P (X)
P (A | c)P (X | c)P (c)
=P
P (c | A, X) = P
c P (A | c)P (X | c)P (c)
c P (A | c)P (c | X)P (X)
P (A | c)P (c | X)
=P
.
c P (A | c)P (c | X)
From (1) and (2), we observe that the distribution P (X) is a “nuisance” in the twostep procedure. Hence we are able to avoid modeling and estimating the marginal
distribution of X. As a result, the effort can be saved for the inference of P (A | c)
and P (c | X).
2.2. Node-coupled Stochastic Block Model. The conditional independence and
follow-up arguments in Section 2.1 pave the way to a flexible framework of models for networks with nodal covariates: specifying the two conditionals P (A | c)
and P (c | X). A similar modeling strategy was proposed in Newman and Clauset
(2015), along with detailed empirical results. Unlike them, we will consider a different model and present a thorough study from both theoretical and computational
perspectives. Note that the conditional distribution P (A | c) only involves the edge
information of the network, while P (c | X) is often encountered in the standard
regression setting for i.i.d. data. This motivates us to consider the following model.
Node-coupled Stochastic Block Model (NSBM):
Q Aij
(a) P (A | c) =
Bci cj (1 − Bci cj )1−Aij ,
i<j

(b) P (c | X) =

Q
i

exp(βcTi xi )
PK
,
T
k=1 exp(βk xi )

where B = (Bab ) ∈ [0, 1]K×K is symmetric, β = (β1 , . . . , βK ) ∈ RKp . The
distribution P (A | c) in (a) follows the stochastic block model (SBM), which, as
the fundamental model, has been extensively studied in the literature. The SBM
implies that the distribution of an edge between node i and j only depends on
their community membership ci and cj . The nodes from the same community are

6

H. WENG AND Y. FENG

stochastically equivalent. The element Bab in the matrix B represents the probability of edge connection between a node in community a and a node in community b.
The P (c | X) in (b) simply takes a multi-logistic regression form, where we will
assume βK = 0 for identifiability. Simple as it looks, we would like to point out
some advantages of NSBM:
• The parameters in NSBM can be estimated by combining the estimation of
B under (a) and β under (b), as we shall elaborate in Section 4.
• The coefficient β reflects the contribution of each nodal covariate for identifying community structures. This information can help us better understand
the implication of the network communities.
• The probability p(c = k | x) =

exp(βkT x)
PK
T
k=1 exp(βk x)

can be used to predict a new

node’s community membership c based on its covariates x, without waiting
for it to form network connections.
• Both P (A | c) and P (c | X) can be readily generalized to fit more complicated structures.
R EMARK 1.

As illustrated in Section 2.1, under the conditional independence

assumption A ⊥ X | c, it is sufficient to consider the conditional likelihood
P (A, c | X) instead of the full version P (A, c, X). In particular, we study the
maximum likelihood estimate, maximum variational likelihood estimate, and the
maximum profile likelihood estimate based on the conditional likelihood in the next
section. However, we emphasize that the conditional independence assumption is
not part of NSBM, though it was used to motivate the model. In the next section, we
will treat this assumption as working independence to derive the likelihood based
estimates. Hence, the three estimates are in fact based on pseudo likelihood. With a
slight abuse of terminology and for simplicity, we still call them the aforementioned
likelihood names in the rest of the paper.
3. Statistical Inference under NSBM. For community detection, our main
goal is to find an accurate community assignment estimator ĉ for the underlying
true communities c. Theoretically, we would like to study the consistency of community detection for a given method. We adopt the notions of consistency from

COMMUNITY DETECTION WITH NODAL INFORMATION

7

Bickel and Chen (2009) and Zhao et al. (2012):
(strong consistency)
(weak consistency)

P (ĉ = c) → 1, as n → ∞,
n
1 X

∀ > 0, P
1(ĉi 6= ci ) <  → 1, as n → ∞.
n
i=1

As the network size increases to infinity, with probability approaching 1, strong
consistency requires perfect recovery of the true community structure, while weak
consistency only needs the mis-classification rate to be arbitrarily small. Note that
since community structure is invariant under a permutation of the community labels
in {1, 2, . . . , K}, the consistency notations above as well as the estimators to be
introduced should always be interpreted up to label permutations.
In the asymptotic setting where the network size n → ∞, holding the parameter B ∈ [0, 1]K×K unchanged implies that the total number of edges present in
the network is of order O(n2 ). Such networks are unrealistically dense. To study
under a more realistic asymptotic framework, we allow B to change with n. In
particular, we consider a sequence of submodels where B = ρn B̄ with B̄ fixed
and ρn = P (Aij = 1) → 0 as n → ∞. The same asymptotic formulation was
studied in Bickel and Chen (2009); Zhao et al. (2012); Bickel et al. (2013). In this
way, the parameter ρn directly represents the sparsity level of the network. For the
consistency results to be derived in the subsequent sections, we will specify the
sufficient conditions on the order of ρn .
As pointed out in Section 2.2, the parameter β in NSBM is associated with the
contribution of each nodal covariate for discovering communities. Measuring the
importance of each nodal attribute to the community structure may provide insightful information about the network. For that purpose, in addition to community
detection, we will study the asymptotics of the estimators for β as well. Since the
parameter B is not of current interest, we will skip the theoretical analysis of the
corresponding estimators.
3.1. Consistency of Maximum Likelihood Method. In Section 2.1, we pointed
out the appealing implication of the assumed conditional independence for likelihood based inference procedure. We now evaluate this procedure under the asymp-

8

H. WENG AND Y. FENG

totic framework we introduced at the beginning of Section 3. Towards that end, we
define the following maximum likelihood based estimators2 :
(3)
(β̂, B̂) =

arg max
βK =0, β∈RKp
B∈[0,1]K×K ,B T =B

X Y Aij
Y exp(βcT xi )
i
,
Bci cj (1 − Bci cj )1−Aij ·
PK
T
exp(β
k=1
c i<j
k xi )
i

(4)
Y Aij
Y exp(β̂cT xi )
i
B̂ci cj (1 − B̂ci cj )1−Aij ·
.
PK
n
exp(
β̂kT xi )
c∈{1,...,K} i<j
k=1
i

ĉ = arg max

The estimators defined above are the realizations of the two-step procedure we
mentioned at the beginning of Section 2.1. We are mainly interested in studying
the consistency of β̂ and ĉ.
First, we would like to introduce several technical conditions.
C ONDITION 1.

B̄ has no two identical columns.

If the probability matrix B̄ has two identical columns, then there exist at least
two communities unidentifiable with each other. In practice, it makes sense to combine those communities into a bigger one.
C ONDITION 2.

iid

(c1 , x1 ), . . . , (cn , xn ) ∼ (c, x) with E(xxT )  0, where  0

represents the matrix being positive definite.
Condition 2 ensures the coefficient vector β is uniquely identifiable.
C ONDITION 3.

There exist constants κ1 and κ2 such that for sufficiently large

t, we have
P (kxk2 > t) ≤ κ1 e−κ2 t .
Condition 3 imposes a sub-exponential tail bound on kxk2 , which is equivalent
to sub-exponential tail assumption on each component of x, via a simple union
2

Recall that the likelihood formulation is the pseudo version as pointed out in Remark 1. Similar

explanations hold for the other two likelihood based methods presented in the subsequent sections.

COMMUNITY DETECTION WITH NODAL INFORMATION

9

bound argument. This covers many different types of covariates like discrete, Gaussian and exponential.
T HEOREM 1.

Assume the data (A, X) follows NSBM and Conditions 1, 2 and

3 hold. In addition, assume

nρn
log n

P (ĉ = c) → 1,

→ ∞ as n → ∞. Then, we have as n → ∞
√

d

n(β̂ − β) → N (0, I −1 (β)),

where I(β) is the Fisher information for the multi-logistic regression problem of
regressing c on X.
The key condition

nρn
log n

→ ∞ requires that the expected degree of every node to

grow faster than the order of log n. The same condition has been used in Bickel and
Chen (2009) and Zhao et al. (2012) to derive strong consistency under SBM. Under
the conditions of the theorem, the maximum likelihood method not only gives us a
strong consistent community assignment estimate ĉ, but also a coefficient estimate
β̂ which is as efficient as if the true label c were known.
3.2. Consistency of Variational Method. The maximum likelihood method studied in Section 3.1 has been shown to have nice theoretical properties. However, the
likelihood function form in (3) renders the computation of the estimators (β̂, B̂)
intractable. In particular, it is computationally infeasible to even evaluate the likelihood function value at a non-degenerate point (when n is not too small), due to the
marginalization over all possible membership assignments. To address this computation issue, we propose a tractable variational method, and demonstrate that it enjoys equally favorable asymptotic properties as the maximum likelihood approach.
This is motivated by the works about variational methods under SBM (Daudin
et al., 2008; Celisse et al., 2012; Bickel et al., 2013). Throughout this section, we
will use the generic symbol P (·) to denote joint distributions and θ = (β, B). To
begin with, recall the well known identity:
log P (A, X; θ) = EQ [log P (A, X, c; θ) − log Q(c)] + D[Q(c) || P (c | A, X; θ)],
where Q(·) denotes any joint distribution of c; the expectation EQ (·) is taken
with respect to c under Q(c); D[·||·] is the Kullback-Leibler divergence. Since

10

H. WENG AND Y. FENG

D[Q(c) || P (c | A, X; θ)] ≥ 0 and the equality holds when Q(c) = P (c |
A, X; θ), it is not hard to verify the following variational equality,
(5)

max log P (A, X; θ) = max EQ [log P (A, X, c; θ) − log Q(c)].
θ

θ,Q(·)

Hence, to compute the maximum likelihood value, we can equivalently solve the
optimization problem on the right hand side of (5). Note that iteratively optimizing
over θ and Q(·) leads to the EM algorithm (Dempster et al., 1977). However, the
calculation of P (c | A, X; θ) at each iteration of EM is computationally intensive.
Instead of optimizing over the full distribution space of Q(·), variational methods
aim to solve an approximate optimization problem, by searching over a subset of
all possible Q(·). In particular, we consider the mean-field variational approach
(Jordan et al., 1999),
max EQ [log P (A, X, c; θ) − log Q(c)],

(6)

θ,Q∈Q

where Q = {Q : Q(c) =

Qn

i=1 qici ,

P

k qik

= 1, 1 ≤ i ≤ n}. The subset Q con-

tains all the distributions under which the elements of c are mutually independent.
The independence structure turns out to make the computation in (6) manageable.
We postpone the detailed calculations to Section 4, and focus on the asymptotic
analysis in this section. Denote the maximizer in (6) by (β̌, B̌) and
(7)

Y exp(β̌cT xi )
Y Aij
i
.
B̌ci cj (1 − B̌ci cj )1−Aij ·
PK
n
exp(
β̌kT xi )
c∈{1,...,K} i<j
k=1
i

č = arg max

T HEOREM 2.

Suppose the conditions in Theorem 1 hold. Then as n → ∞
P (č = c) → 1,

√

d

n(β̌ − β) → N (0, I −1 (β)),

where I(β) is the Fisher information for the multi-logistic regression problem of
regressing c on X.
As we can see, under the same conditions as the maximum likelihood method,
the variational approach can deliver equally good estimators, at least in the asymptotic sense. In other words, the approximation made by the variational method does

COMMUNITY DETECTION WITH NODAL INFORMATION

11

not degrade the asymptotic performance. This should be attributed to the condition
nρn
log n

→ ∞, which guarantees the network has sufficient edge information for doing

approximate inference.
3.3. Consistency of Maximum Profile Likelihood Method. The two methods
presented in Sections 3.1 and 3.2 are implementations of the two-step procedure
we discussed in Section 2.1: first estimating parameters based on the likelihood
function and then doing posterior inference using the estimated distribution. In this
section, we introduce a one-step method that outputs the parameter and community
assignment estimates simultaneously. The method solves the following problem,
(8)
(β̃, B̃, c̃) =

arg max
βK =0, β∈RKp
c∈{1,...,K}n
B∈[0,1]K×K ,B T =B

Y exp(βcT xi )
Y Aij
i
.
Bci cj (1 − Bci cj )1−Aij ·
PK
T
exp(β
k=1
k xi )
i
i<j

In the above formulation, we treat the latent variables c as parameters and obtain
the estimators as the maximizer of the joint likelihood function. This enables us
to avoid the cumbersome marginalization encountered in the maximum likelihood
method. This approach is known as maximum profile likelihood (Bickel and Chen,
2009; Zhao et al., 2012). Bickel and Chen (2009) showed strong consistency under
stochastic block model, and Zhao et al. (2012) generalized the results to degreecorrected block models. Following similar ideas, we will investigate this method in
the node-coupled stochastic block model. For theoretical convenience, we consider
a slightly different formulation:
(9)

(β̃, B̃, c̃) =

arg max
βK =0, β∈RKp
c∈{1,...,K}n
B∈RK×K ,B T =B

Y −B
Y exp(βcT xi )
A
i
,
e ci cj Bci cijj ·
PK
Tx )
exp(β
i
k=1
k
i<j
i

where the Bernoulli distribution in (8) is replaced by Poisson distribution. In our
asymptotic setting ρn → 0, the difference becomes negligible.
T HEOREM 3.
hold.

Assume the data (A, X) follows NSBM and Conditions 1 and 2

12

H. WENG AND Y. FENG

(i) If nρn → ∞ and Ekxkα2 < ∞ (α > 1), then there exists a constant γ > 0
such that, as n → ∞
P

n
1 X

n



1(c̃i 6= ci ) ≤ γ(nρn )−1/2 → 1, kβ̃ − βk2 = Op ((nρn )

1−α
2α

).

i=1

(ii) Assume Condition 3 is satisfied. If
P (c̃ = c) → 1,

√

nρn
log n

→ ∞, then as n → ∞
d

n(β̃ − β) → N (0, I −1 (β)),

where I(β) is the Fisher information for the multi-logistic regression problem of regressing c on X.
We see that part (ii) in Theorem 3 is identical to Theorems 1 and 2. Hence the
maximum profile likelihood method is equivalently good as the previous two, in
certain sense. The conclusions in part (i) shed lights on how the network edges and
nodal covariates affect the consistency results. Under the scaling nρn → ∞, c̃ is
only weak consistent. And the higher moment kxk2 has, the faster convergence rate
√
β̌ can achieve. Suppose all moments of kxk2 exist, then we would have nρn kβ̌−
βk2 = Op (1). Since ρn → 0, this convergence rate is slower than and may be
arbitrarily close to the one in part (ii) when

nρn
log n

→ ∞.

4. Practical Algorithms. In Section 3, we have studied three likelihood based
community detection methods and shown their superb asymptotic performances. In
this section, we design and analyze specialized algorithms, for computing the variational estimators defined by (6), (7) and the maximum profile likelihood estimators
in (9). As discussed in Section 3.2, the maximum likelihood estimators are computationally infeasible, hence omitted here. The key challenge lies on the fact that
the likelihood based functions in (6), (7) and (9) are all non-convex. Multiple local
optima may exist and the global solution is often impossible to accurately allocate.
To address this issue, we first obtain a “well behaved” preliminary estimator via
convex optimization, and then feed it as an initialization into “coordinate” ascent
type iterative schemes. The idea is that the carefully chosen initialization may help
the followed-up iterations to escape “bad” local optima and arrive “closer” (better

COMMUNITY DETECTION WITH NODAL INFORMATION

13

approximation) to the ideal global solution. As we shall see in the numerical studies, the results with a “well behaved” initial estimator are significantly better than
those with a random initialization. These two steps will be discussed in detail in
Sections 4.1 and 4.2, respectively.
4.1. Initialization via Convex Optimization. The convex optimization we consider in this section is semidefinite programming (SDP). Different formulations of
SDP have been shown to yield good community detection performances in Chen
et al. (2012); Amini and Levina (2014); Cai and Li (2015); Montanari and Sen
(2015); Guédon and Vershynin (2016), among others. One illuminating interpretation of SDP is to think of it as a convex relaxation of the maximum likelihood
method. For example, starting from a specialized stochastic block model, one can
derive SDP by approximating the corresponding likelihood function. See Chen
et al. (2012); Amini and Levina (2014); Cai and Li (2015) for the detailed arguments. However, under NSBM, because of the nodal covariates term, it is not
straightforward to generalize the convex relaxation arguments. We hence resort to
a different understanding of SDP elaborated in Guédon and Vershynin (2016). The
key idea is to construct SDP based on the observations directly, with the goal of
having the true community assignment c to be the solution of a “population” version of the SDP under construction. Then with a few conditions, we would like to
show that the solution of the SDP is “close” to the solution of its “population” version, i.e., the true community assignment. In particular, we consider the following
semidefinite programming problem,
(10)

Ẑ = arg maxhA + γn XX T , Zi
Z

subject to Z  0, Z ∈ Rn×n
0 ≤ Zij ≤ 1, 1 ≤ i, j ≤ n
X
Zij = λn ,
ij

where γn , λn > 0 are two tuning parameters; h·, ·i denotes the inner product of
two matrices. We then obtain the communities by running K-means on Ẑ (treating

14

H. WENG AND Y. FENG

each row of Ẑ as a data point in Rn ). We show that the approach can produce a
consistent community assignment estimate as presented in the following theorem.
T HEOREM 4.

iid

Assume part (a) in NSBM holds and (c1 , x1 ), . . . , (cn , xn ) ∼

(c, x). Let {c̄i }ni=1 be the community estimates from running K-means on Ẑ defined
in (10). If mina B̄aa > maxa6=b B̄ab , nρn → ∞ and kxk2 is sub-Gaussian, then by
Pn
P
2
choosing γn = o(ρn ) and λn = K
k=1 ( i=1 1(ci = k)) , we have
n

1X
P
1(c̄i 6= ci ) →
0.
n
i=1

The proof of Theorem 4 will provide a clear picture on why (10) is constructed
in that way. But to avoid digression, we defer the proof to the Appendix. As seen
from Theorem 4, though the SDP works under more assumptions than previously
discussed likelihood based methods, the conditions are not very stringent. The crucial assumption mina B̄aa > maxa6=b B̄ab requires denser edge connections within
communities than between them, which is satisfied by most real networks. Furthermore, the two tuning parameters (γn , λn ) coupled with (10) need to be chosen
appropriately. The tuning γn trades off the information from two different sources:
network edge and nodal covariates. The way we incorporate nodal covariates has
the same spirit as Binkiewicz et al. (2014) does in spectral clustering. From the
simulation studies in the next section, we shall see that a flexible choice of γn can
lead to satisfactory results. The choice of parameter λn in Theorem 4 depends on
the unknown truth in a seemingly restrictive way. However, we will demonstrate
through simulations that the community detection results are quite robust to the
choice of λn .
The convex optimization problem (10) can be readily solved by standard semidefinite programming solvers such as SDPT3 (Tütüncü et al., 2003). However, those
solvers are based on interior-point methods, and are computationally expensive
when the network size n is more than a few hundred. To overcome this limit, we
apply the alternating direction method of multipliers (ADMM) to develop a more
scalable algorithm for solving (10). We start by a brief description of the generic
ADMM algorithm with the details available in the excellent tutorial by Boyd et al.

15

COMMUNITY DETECTION WITH NODAL INFORMATION

(2011). In general, ADMM solves problems in the form
minimize f (y) + h(z)

(11)

subject to By + Dz = w,
where y ∈ Rm , z ∈ Rn , B ∈ Rq×m , D ∈ Rq×n , w ∈ Rq ; and f (y), h(z) are two
convex functions. The algorithm takes the following iterations at step t.
(12)
(13)



y t+1 = arg min f (y) + (ξ/2)kBy + Dz t − w + ut k22 ,
y


z t+1 = arg min h(z) + (ξ/2)kBy t+1 + Dz − w + ut k22 ,
z

(14)

t+1

u

t

= u + By t+1 + Dz t+1 − w,

with ξ > 0 being a step size constant.
To use this framework, we reformulate (10) as:
minimize l(Z  0) + l(0 ≤ Yij ≤ 1, 1 ≤ i, j ≤ n) + l

X

Wij = λn



ij

− hA + γn XX T , Zi
subject to Y = Z, Y = W,
where Z, Y, W ∈ Rn×n ; l(Z  0) equals 0 if Z  0 and +∞ otherwise; similar
2

definitions hold for other l(·). If we set y = (vec(Y ), vec(Y ))T ∈ R2n , z =
2

2 ×2n2

(vec(W ), vec(Z))T ∈ R2n , B = −D = I2n2 ∈ R2n

2

, w = 0 ∈ R2n , where

vec(·) denotes the vectorized version of a matrix, then the problem above becomes
an instance of (11). The corresponding iterations have the following expressions:


Y t+1 = arg min kY − W t + U t k2F + kY − Z t + V t k2F ,
0≤Yij ≤1

W t+1 = Parg min kY t+1 − W + U t k2F ,
ij

Wij =λn



Z t+1 = arg min − hA + γn XX T , Zi + (ξ/2)kY t+1 − Z + V t k2F ,
Z0

U t+1 = U t + Y t+1 − W t+1 , V t+1 = V t + Y t+1 − Z t+1 ,

16

H. WENG AND Y. FENG

where k · kF denotes the Frobenius norm. It is not hard to see that each iteration
above has a closed form update with the details summarized in Algorithm 1.3
Algorithm 1 Solving (10) via ADMM
Input: initialize Z 0 = A + γn XX T , W 0 = Y 0 = U 0 = V 0 = 0, number of iterations T , step
size ξ.
For t = 0, . . . , T − 1
(a) Y t+1 = min{max{0, 12 (W t + Z t − U t − V t )}, 1}.
P
t
(b) W t+1 = Y t+1 + U t + n−2 [λn − ij (Yijt+1 + Uij
)]11T .
(c) Z t+1 = P Λ+ P T , where Y t+1 + V t + ξ −1 · (A + γn XX T ) = P ΛP T .
(d) U t+1 = U t + Y t+1 − W t+1 , V t+1 = V t + Y t+1 − Z t+1 .
Output Z T .

4.2. Coordinate Ascent Scheme. As we may see, the problem formulations in
(6) and (9) are not suitable for gradient or Hessian based iterative algorithms, because they either involve discrete variables or have non-trivial constraints. The variables involved in those optimization problems can be divided into community
assignment related and others. Naturally, we will adopt the iterative scheme that
alternates between these two types of variables.
4.2.1. Computing Variational Estimates. To compute the variational estimates
in (6), we follow the EM style iterative fashion by maximizing the objective function in (6) with respect to θ and Q ∈ Q alternatively. Specifically, we are solving

(15)

Xh X

β t+1 = arg max
β∈RpK ,βK =0

i

t
qik
βk

k

T

xi − log

X

T

eβk xi

i
,

k

(16)
B t+1 =

arg max

Xh

B∈RK×K ,B T =B ab
3

log Bab ·

X
i<j

t t
Aij qia
qjb + log(1 − Bab ) ·

X

i
t t
(1 − Aij )qia
qjb ,

i<j

In Step (c), P ΛP T denotes the spectral decomposition; Λ+ represents the truncated (keep posi-

tive elements) version of Λ.

17

COMMUNITY DETECTION WITH NODAL INFORMATION

(17)
t+1
{qik
} = arg max
{qik }

Xh

t+1
log Bab
·

X

t+1
Aij qia qjb + log(1 − Bab
)·

i<j

ab

X

(1

i<j

i XX
XX
− Aij )qia qjb +
qik (βkt+1 )T xi −
qik log qik .
i

i

k

k

Note that the objective function in (15) takes a similar form as the log-likelihood
function of multi-logistic regression model. We hence use Newton-Raphson algorithm in the same way as we fit multi-logistic regression model, to compute the
update in (15). This corresponds to Step (a) of Algorithm 2, in which we have used
the name “FitMultiLogistic” there to denote the full step with a bit abuse of notation. In addition, the update in (16) has an explicit solution, which corresponds to
Step (b) in Algorithm 2. Regarding the update for {qik }ik in (17), unfortunately,
the optimization is non-convex and does not have analytical solutions. We then
implement an inner blockwise coordinate ascent loop to solve it. In particular, we
update {qik }K
k=1 one at a time:
hXX
X
qik ·
Aij qjb · log Bkb + (1 − Aij )qjb · log(1
{qik }k = arg max
{qik }k

k

b

− Bkb )

i

+

j6=i

X

qik βkT xi −

k

X

qik log qik .

k

It is straightforward to show that the update above has closed forms:


XX
eak
, ak = βkT xi +
qjb · Aij log Bkb + (1 − Aij ) log(1 − Bkb ) .
qik = PK
ak
k=1 e
b j6=i
T }, β T , B T via AlgoThis yields Step (c) for Algorithm 2. After computing {qik

rithm 2, we calculate the community assignment estimate č based on (7). This
could be done by coordinate ascent iterations, like Step (c) in Algorithm 3 (to be
introduced in Section 4.2.2. Alternatively, we can use the following approximated
T }:
posterior distribution {qik

(18)

T
či = arg max qik
,

1 ≤ i ≤ n.

1≤k≤K

In numerical studies, we adopt the approach in (18), which is computationally more
efficient.

18

H. WENG AND Y. FENG

Algorithm 2 Solving (6) via iterating between (β, B) and Q.
0
Input: initialize {qik
}, number of iterations T
For t = 0, . . . , T − 1
t
(a) β t+1 = FitMultiLogistic(X, {qik
})
t+1
(b) Bab
=

P

i<j
P

t
t
Aij qia
qjb

i<j

t qt
qia
jb

t+1
(c) Update {qik
} via Repeating

For i = 1, . . . , n

P P
t+1
t+1
log qik ∝ (βkt+1 )T xi + b j6=i qjb · Aij log Bkb
+(1−Aij ) log(1−Bkb
)
T
Output {qik
}, β T , B T .

4.2.2. Computing Maximum Profile Likelihood Estimates. Similarly to the variational estimates, we maximize the likelihood function in (9) with respect to (β, B)
and c iteratively. In other words, we solve
β

(19)

t+1

Xh
X T i
T
= arg max
βct xi − log
eβk xi ,
β∈RpK ,βK =0

i

i

k

(20)
t+1
Bab
= arg max log Bab ·
Bab

(21)

X

Aij 1(cti = a, ctj = b) − Bab ·

i<j

ct+1 = arg max

Xh
X
i<j

1(cti = a, ctj = b),

i<j

c∈{1,...,K}n ab
t+1
·
− Bab

X

t+1
log Bab
·

X

Aij 1(ci = a, cj = b)

i<j

i

1(ci = a, cj = b) +

X
(βct+1
)T xi .
i
i

Here, solving (19) is equivalent to computing the maximum likelihood estimate of
multi-logistic regression. This is carried out in Step (a) of Algorithm 3. In addition,
it is straightforward to see that Step (b) in Algorithm 3 is the solution to (20). For
computing ct+1 in (21), we update its element one by one, as shown by Step (c) in
Algorithm 3.
4.2.3. Variational Estimates vs. Maximum Profile Likelihood Estimates. So far
we have studied the theoretical properties of variational and maximum profile like-

COMMUNITY DETECTION WITH NODAL INFORMATION

19

Algorithm 3 Solving (9) via iterating between (β, B) and c.
Input: initialize {c0i }, number of iterations T
For t = 0, . . . , T − 1
(a) β t+1 = FitMultiLogistic(X, ct )
P
t
t
t+1
i<j Aij 1(ci =a,cj =b)
(b) Bab
= P
t
t
i<j 1(ci =a,cj =b)
(c) Update {ct+1
} via Repeating
i
For i = 1, . . . , n
P P
t+1
t+1
ci = arg max1≤k≤K (βkt+1 )T xi + b j6=i 1(cj = b)·(Aij log Bkb
−Bkb
)
Output {cTi }, β T .

lihood estimates, and developed algorithms to compute them. The results in Sections 3.2 and 3.3 demonstrate that they have the same asymptotic performance under

nρn
log n

→ ∞. We now compare the corresponding algorithms. By taking a close

look at Algorithms 2 and 3, we observe that the three steps in the two algorithms
share a lot of similarities. Algorithm 2 is essentially a “soft” version of Algorithm 3
in the following sense: instead of using the community assignment ci in Algorithm
3, the steps in Algorithm 2 involve the probability of belonging to every possible
community. This might remind us of the comparison between the EM algorithm
and K-means under Gaussian mixture models. As we will see in Section 5, variational and maximum profile likelihood methods usually lead to similar numerical
results.
5. Numerical Experiments. In this section, we conduct a detailed experimental study of the SDP defined in (10), variational and maximum profile likelihood
methods on both simulated and real datasets. We use two quantitative measures for
evaluating their community detection performance.
Normalized Mutual Information (Ana and Jain, 2003):


P P
n ·n
−2 i j nij log ni·ijn·j
  P
 .
NMI = P
n·j
ni·
n
log
+
n
log
i·
·j
i
j
n
n

20

H. WENG AND Y. FENG

Adjusted Rand Index (Hubert and Arabie, 1985):
(n2i· ) j (n2·j )
−
(n2 )
ARI = P
 1 P n·j  Pi (n2i· ) Pj (n·j ) .
ni·
1
2
i 2 + 2
j 2 −
2
(n2 )
nij 
ij 2

P

P

i

P

In the above expressions, ni· denotes the true number of nodes in community i, n·j
represents the number of nodes in the estimated community j and nij is the number
of nodes belonging to community i but estimated to be in community j. Both NMI
and ARI are bounded by 1, with the value of 1 indicating perfect recovery while
0 implying the estimation is no better than random guess. See Steinhaeuser and
Chawla (2010) for a detailed discussion.
1.5

5.1. Simulation Studies. We set K = 2, ρn = 3[log(n)]
, P (c = 1) = P (c =
4n
!
1.6 0.4
. We consider the following two different scenarios.
2) = 0.5, B̄ =
0.4 1.6
(A). p = 4, x | c = 1 ∼ N (µ, I4 ), x | c = 2 ∼ N (−µ, I4 ), µ = (0, 0.4, 0.6, 0.8)T ,
where I4 ∈ R4×4 is the identity matrix.
(B). p = 4, (x1 , x2 ) | c = 1 ∼ N (µ, Σ), (x1 , x2 ) | c = 2 ∼ N (−µ, Σ), µ =
(0.5, 0.5)T , Σ11 = Σ22 = 1, Σ12 = 0.3, x3 | c = 1 ∼ Bernoulli(0.6), x3 |
c = 2 ∼ Bernoulli(0.4), x4 | c = 1 ∼ Uniform(−0.2, 0.5), x4 | c = 2 ∼
Uniform(−0.5, 0.2); and (x1 , x2 ), x3 and x4 are mutually independent.
Note that in Scenario (A), NSBM is the correct model and the first nodal variable is independent of the community assignment; In Scenario (B), NSBM is no
longer correct. Under both correct and misspecified models, we would like to: (i)
investigate the impacts of the two tuning parameters (γn , λn ) in the SDP (10); (ii)
examine the effectiveness of the SDP as initialization; (iii) check the performances
of variational and maximum profile likelihood methods for utilizing nodal information.
5.1.1. Tuning Parameters in SDP. For both simulation settings, we solve SDP
defined in (10) with different tuning parameters via Algorithm 1, with the number
of iterations T = 100 and the step size ξ = 1. We then calculate the NMI of

21

COMMUNITY DETECTION WITH NODAL INFORMATION

its community detection estimates. Since the ARI gives similar results, we do not
show them here for simplicity. The full procedure is repeated 500 times.
Figure 1 demonstrates the joint impact of the tuning parameters on the SDP performance under Scenario (A). First of all, the comparison of NMI between γn = 0
and γn > 0 indicates the effectiveness of SDP (10) for leveraging nodal information. We can also see that neither small or large values of γn lead to optimal performances, verifying the point we discussed in Section 4.1 that γn plays the role of
balancing the edge and nodal information. An appropriate choice, as suggested by
the four plots, is γn =

[log(n)]0.5
,
n

which is consistent with the result of Theorem 4.

Regarding the parameter λn , we know from Theorem 4 that λn =

n2
2

is the desired

choice. Interestingly, Figure 1 shows that a wide range of λn can give competitive
results, as long as the corresponding γn is properly chosen. For Scenario (B), similar phenomena can be observed in Figure 2. Note that since the nodal covariates are
not as informative as in Scenario (A), the optimal γn ≈

0.8[log(n)]0.5
n

tends to give

more weights to the adjacency matrix. The results in these two different settings
confirm the implication of Theorem 4, that SDP (10) can work beyond NSBM.
5.1.2. Community Detection Performance via Variational and Maximum Profile
Likelihood Methods. We implement variational and maximum profile likelihood
methods via Algorithms 2 and 3 respectively, taking the outputs from Algorithm 1
as initialization (called VEM-C and MPL-C respectively). We do not predefine the
number of iterations T in both algorithms, and instead keep iterating until convergence. To investigate the impact of SDP as an initialization, we have additionally
implemented both methods with random initialization (called VEM-B and MPLB respectively): run Algorithms 2 and 3 with random initialization independently
multiple times and choose the outputs that give the largest objective function value
(e.g., the profile likelihood function). We have also applied both methods for the
simulated datasets with nodal attributes removed (called VEM-A and MPL-A respectively). This will be used as a comparison to check the effect of the two meth[log(n)]0.5
for all
n
0.8[log(n)]0.5
1 2
2 n , γn =
n

ods in incorporating nodal information. We set λn = 21 n2 , γn =
the implementations of SDP under Scenario (A); and λn =

22

H. WENG AND Y. FENG

0.7
0.63

0.88
0.63

0.80

0.56

0.75

0.49

0.49

0.80
0.78

0.42

0.65

0.36

0.84
0.82

0.70

0.42

0.86

0.56

τ

τ

0.7

0.85

0.76
0.36

0.60

0.74
0.72

0.6

1.3

2.1

2.8

3.6

4.3

0.3

α (n=100)

0.7

1.1

1.5

1.9

2.3

α (n=200)

0.7

0.7
0.90

0.56

0.86

0.49

0.84

0.42

0.82
0.80

0.36

0.90

0.63

0.88

τ

τ

0.63

0.56

0.88

0.49

0.86

0.42

0.84

0.36

0.82

0.78
0.2

0.5

0.8

1

α (n=300)

1.3

1.6

0.2

0.4

0.6

0.8

1

1.2

α (n=400)

F IG 1. The community detection performance of SDP (measured by NMI), under Scenario (A), with
different tuning parameters (λn , γn ); NMI is averaged over 500 repetitions; We have used the scaled
version of the tuning parameters: τ = λnn2 , α = 100γn .

under Scenario (B).
Figure 3 shows the community detection results of both methods under Scenario
(A). By comparing the four curves in each plot, we can make a list of interesting
observations: (1) SDP is a good initialization (MPL-C vs. MPL-B, VEM-C vs.
VEM-B); (2) SDP itself already gives reasonable outputs, but the follow-up iterations further improve the performance (SDP vs.. MPL-C, SDP vs. VEM-C); (3)
the nodal covariates are helpful for detecting communities, and the two methods
have made effective use of it (MPL-A vs. MPL-C, VEM-A vs. VEM-C); (4) the
two methods have similar performances when initialized with SDP (MPL-A vs.
VEM-A, MPL-C vs. VEM-C). Moreover, we would like to point out the different
behavior of the two methods with random initialization. The comparison between
the two purple curves (MPL-B vs. VEM-B) implies that compared to the variational method, maximum profile likelihood method has the potential of exploring

23

COMMUNITY DETECTION WITH NODAL INFORMATION

0.7

0.7
0.86

0.80

0.63

0.63

0.75

0.56

0.84
0.82

0.56

0.70

0.42

0.65

0.42

0.36

0.60

0.36

τ

τ

0.80
0.49

0.49

0.78
0.76
0.74

0.6

1.2

1.9

2.6

3.3

4

0.72
0.3

α (n=100)

0.7

1.1

1.5

1.9

2.3

α (n=200)

0.7

0.7

0.90

0.88
0.63

0.86
0.84
0.82

0.49

0.86

0.56

0.84

τ

τ

0.56

0.88

0.63

0.49
0.82

0.80
0.42

0.78
0.76

0.36

0.42

0.80

0.36

0.78

0.74
0.2

0.5

0.8

1

1.3

1.6

0.76
0.2

α (n=300)

0.4

0.6

0.8

1

1.2

α (n=400)

F IG 2. The community detection performance of SDP (measured by NMI), under Scenario (B), with
different tuning parameters (λn , γn ); NMI is averaged over 500 repetitions; We have used the scaled
version of the tuning parameters: τ = λnn2 , α = 100γn .

the parameter space more efficiently, especially when the sample size is large. One
possible explanation is that the update of the “soft” community labels (the distribution {qik }) in the variational algorithm may cause it to move very slowly in the
parameter space and hence it may take many steps to change a label assignment.
Furthermore, note that we can use the asymptotic normality property of the estimators for β in Theorems 2 and 3 to perform variable selection. The results of the
Wald test regarding each component of β are presented in Figure 4. We see that
our methods are able to identify relevant (the last three) and irrelevant (the first)
nodal variables. Regarding Scenario (B), similar observations on the community
detection performance can be made from Figure 6. We thus omit the details. As
a final remark, the performances in Scenario (B) indicate that both methods can
work to a certain extent of model misspecification.

24

H. WENG AND Y. FENG

NMI

Variational EM

SDP
MPL-A
MPL-B
MPL-C

100

200

300

sample size

400

500

0.70 0.75 0.80 0.85 0.90

NMI

0.70 0.75 0.80 0.85 0.90

Maximum Profile Likelihood

SDP
VEM-A
VEM-B
VEM-C

100

200

300

400

sample size

500

F IG 3. The community detection performances under Scenario (A); The average NMI is calculated
over 500 repetitions along with its standard error bar; MPL-A, MPL-B and MPL-C represent the
maximum profile likelihood methods with no nodal covariates being used, random initialization, and
initialization from SDP, respectively; Similar notations are used for variational method. We have
used 15 independent random initializations for maximum profile likelihood method across all the
sample size; For variational method, the number of random initializations used starts from 15 for
n = 100 and consecutively increases by 1 for the subsequent sample sizes.

5.2. Real Data Analysis. The dataset is about a research team consisting of
77 employees in a manufacturing company (Cross and Parker, 2004). A weight
wij is assigned to the edge from employee i to employee j based on the extent
to which employee i provides employee j with information j uses to accomplish
j’s work. There are seven choices for the weights: 0 (I do not know this person/I
have never met this person); 1 (Very infrequently); 2 (Infrequently); 3 (Somewhat
infrequently); 4 (Somewhat frequently); 5 (Frequently); 6 (Very frequently). In addition to the edge information, the dataset also contains several attributes of each
employee: location (1: Paris, 2: Frankfurt, 3: Warsaw, 4: Geneva); tenure (1: 112 months, 2: 13-36 months, 3: 37-60 months, 4: 61+ months); the organizational
level (1: Global Dept Manager, 2: Local Dept Manager, 3: Project Leader, 4: Researcher). Since the network is a weighted and directed network, we first convert
it to a binary network such that there exists an edge from i to j if and only if

25

1.0

COMMUNITY DETECTION WITH NODAL INFORMATION

0.8
0.6

Power

0.04
0.02

β2
β3
β4

0.4

0.00

Type I error

β1

100

200

300

sample size

400

500

100

200

300

sample size

400

500

F IG 4. Wald test for each component of β. The calculations are averaged over 2000 repetitions. The
significance level is set to 0.01. Since both variational and maximum profile likelihood methods give
similar results, we only present the result of variational method for simplicity.

wij > 3. This corresponds to whether the information is provided frequently or
not. We then further convert it into an undirected network in the way that the edge
between i and j exists if and only if both directed edges from i to j and j to
i are present. Finally, we remove three isolated nodes from the network. To explore the intro-organizational community structure, we re-order the adjacency matrix based on random permutation and the attributes. As can be seen from Figure 5,
the attribute “location” is a very informative indicator of the network’s community
structure. This should not come as a big surprise, since the same office location

random

location

tenure

level

F IG 5. From left to right are the re-ordered adjacency matrices based on random permutation, location, tenure, and organizational level.

26

H. WENG AND Y. FENG

NMI

Variational EM

SDP
MPL-A
MPL-B
MPL-C

100

200

300

400

sample size

500

0.70 0.75 0.80 0.85 0.90

NMI

0.70 0.75 0.80 0.85 0.90

Maximum Profile Likelihood

SDP
VEM-A
VEM-B
VEM-C

100

200

300

400

sample size

500

F IG 6. The community detection performances under Scenario (B); All the relevant descriptions are
the same as in Figure 3.

usually promotes communication and collaboration between team members. We
now use the “location” as the ground truth for the community assignment and examine the performances of SDP (10), maximum profile likelihood and variational
methods based on the rest of the data. For SDP (10), we first use spectral clustering on adjacency matrix (Lei and Rinaldo, 2014) to estimate the size of the communities and plug the estimates in the formula of λn in Theorem 4 to determine
ρ̂n
λn . Regarding γn , motivated from the simulation results, we choose γn = log
n,
2×number of edges
where ρ̂n =
. The maximum profile likelihood and variational
n2

methods are initialized by the output from SDP. We can see from Table 1 that by
incorporating the nodal information, community detection accuracy has been improved. It is interesting to observe that SDP performs as well as the two likelihood
based methods, when nodal covariates are available. Note that we can calculate the
mutual information between the “ground truth” variable “location” and the other
two to see how much community information they contain. Given that both mutual
information (0.11 & 0.03) are pretty small, the magnitude of improvement in Table
1 is reasonable.

COMMUNITY DETECTION WITH NODAL INFORMATION

edge
edge + nodal

SDP
0.881
0.920

MPL
0.894
0.920

VEM
0.894
0.920

SDP
0.882
0.921

MPL
0.883
0.921

27

VEM
0.883
0.921

TABLE 1
The community detection results of SDP, maximum profile likelihood, and variational methods.
MPL and VEM denote maximum profile likelihood and variational methods, respectively. NMI is
computed on the left part of the table, and ARI on the right. The row indexed by “edge” shows the
results based on the network without nodal information, while the other one “edge+nodal”
contains the results of making use of the two attributes available.

6. Discussion. In this paper, we present a systematic study of the community
detection with nodal information problem. We propose a flexible network modeling framework, and analyze three likelihood based methods under a specialized
model. Both asymptotic and algorithmic aspects have been thoroughly discussed.
The superiority of variational and maximum profile likelihood methods are verified
through a variety of numerical experiments. Finally, we would like to highlight several potential extensions and open problems for future work.
1. The modeling of both the network and nodal covariates can be readily extended to more general families, such as degree-corrected stochastic block
model and non-parametric regression, respectively. The corresponding asymptotic results might be derived accordingly.
2. In the setting with high dimensional covariates, penalized likelihood methods are more appealing for both community detection and variable selection.
Theoretical analysis of community detection and variable selection consistency will be necessary.
3. For very sparse networks, considering nρn = O(1) seems to be a more
realistic asymptotic framework. Under such asymptotic setting, community
detection consistency is impossible. The effect of nodal covariates becomes
more critical. It is of great interest to characterize the impact of the nodal
information on community detection.
4. In this work, we assume the number of communities K is known. How to select K is an important problem in community detection. Some recent efforts
towards this direction include Saldana et al. (2015); Le and Levina (2015);

28

H. WENG AND Y. FENG

Wang and Bickel (2015); Lei (2016).
APPENDIX A: APPENDIX
Notations and Preliminaries. Before the proofs, we first introduce some necessary
P
notations. Let πa = P (c = a), π̂a = n1 ni=1 1(ci = a), 1 ≤ a ≤ K. Given a community assignment e ∈ {1, . . . , K}n , define O(e), V (e), T (e), T̂ (e) ∈ RK×K
and f 0 (e), fˆ(e) ∈ RK with their elements being
P
X
1(ei = a, ci = b)
,
Oab (e) =
Aij 1(ei = a, ej = b), Vab (e) = i P
i 1(ci = b)
ij
X
X
π̂a π̂b B̄ab Vka (e)Vlb (e),
πa πb B̄ab Vka (e)Vlb (e), T̂kl (e) =
Tkl (e) =
ab

ab

fk0 (e) =

X
a

Vka (e)πa ,

fˆk (e) =

X
a

Vka (e)π̂a =

nk (e)
.
n

Here, Oab (e) represents the number of edges between communities a and b under
assignment e (twice the number when a = b), Vab (e) represents the proportion
of the nodes in community b under c that are mislabeled as community a under
e, Tkl (e) reflects the connection probability (up to a scaling factor ρn ) between a
node in community k and another node in community l under e with T̂kl (e) being
its empirical version, and fˆk (e) represents the proportion of community k under e
P
with fk0 (e) as its “population” version. Also denote na (e) = i 1(ei = a), 1 ≤


P 
P
;
H(T,
f
)
=
−
T
a ≤ K; F (T, f ) = ab Tab log fTaab
ab
ab Tab log B̄ab −
fb

P
fa fb B̄ab ; V = {V ∈ RK×K : k Vka = 1, Vka ≥ 0, 1 ≤ k, a ≤ K}; µn =
n2 ρn . Throughout the proofs, we use B 0 and β0 to represent the true parameters
in NSBM while leaving B and β as the generic parameter. We will C1 , C2 , . . . , to
represent positive generic constants whose values may vary across different lines.
P
We will frequently use the notation kV (e) − V (c)k1 = ab |Vab (e) − Vab (c)|.
P
Note that n1 i 1(ei 6= ci ) ≤ 12 kV (e) − V (c)k1 (see the derivation on Page 22
of Zhao et al. (2012)). Moreover, we cite two concentration inequality results that
will be used in the proof. The first one is Lemma A.1 from Zhao et al. (2012):
(22)



Oab (e)
2 µn
P max max
− T̂ab (e) ≥  ≤ 2K n+2 exp(−
),
e
ab
µn
8 maxab B̄ab

COMMUNITY DETECTION WITH NODAL INFORMATION

29

for  < 3 maxab B̄ab . The second one is (1.4) in Bickel et al. (2015): ∀ γ0 > 0,
(23)

P

max

0<kV (e)−V (c)k1 ≤δn

where Xab (e) =


[max |Xab (e) − Xab (c)| − γ0 · kV (e) − V (c)k1 ] ≤ 0 → 1,
ab

Oab (e)
µn −T̂ab (e) and δn

→ 0, nρn / log n → ∞. Finally, we should

be aware that the expressions involving community assignment e in the proofs
are to be interpreted, up to permutations of community labels in {1, 2, . . . , K}
whenever necessary.

L EMMA 1.

There exist positive constants c0 , c1 , c2 , c3 > 0, such that

F (T (c), f0 (c)) − F (T (e), f0 (e)) ≥ c0 · kV (e) − V (c)k1 , if kV (e) − V (c)k1 ≤ c1 ,
H(T (c), f0 (c)) − H(T (e), f0 (e)) ≥ c2 · kV (e) − V (c)k1 , if kV (e) − V (c)k1 ≤ c3 .
P ROOF. We only show the proof for the first inequality since the second one
can be derived in a similar way. Note that F (T (e), f0 (e)) can be considered as
a function of V (e). We give it another notation H(V (e)), where H(·) is defined
on the convex set V. Further define g(; V (e)) = H((1 − )V (c) + V (e)), for
0 ≤  ≤ 1. Since V (e), V (c) ∈ V, g(; V (e)) is well defined. We first show that
∃ c̃0 , c̃1 > 0, s.t.
(24)

g 0 (; V (e)) ≤ −c̃0 · kV (e) − V (c)k1 , for any 0 ≤  ≤ c̃1 , V (e) ∈ V

where g 0 (; V (e)) is the derivative with respect to ; the constants c̃0 , c̃1 do not
depend on V (e). To prove (24), denote Ṽka = (1 − )Vka (c) + Vka (e), rkl =
P
ab πa πb Ṽka Ṽlb B̄ab
P
.
ab πa πb Ṽka Ṽlb

Then a straightforward calculation gives us,

(25)
g 0 (; V (e)) = 2

XX
kl

+

XX
kl

ab

πa πb · [Vka (e) − Vka (c)] · [Vlb (e) − Vlb (c)] · (B̄ab log rkl − rkl )

ab

πa πb ([Vka (e) − Vka (c)]Vlb (c) + [Vlb (e) − Vlb (c)]Vka (c)) · (B̄ab log rkl − rkl ).

30

H. WENG AND Y. FENG

Hence when  = 0, the above equation can be simplified as
(26)
g 0 (0; V (e)) =

XX
ka

+

XX
lb

=

ka

πa πb · (Vka (e) − Vka (c)) · (B̄ab log B̄kb − B̄kb )

b

XX
lb

=2

πa πb · [Vlb (e) − Vlb (c)] · Vka (c) · (B̄ab log B̄kl − B̄kl )

ka

XX
+

πa πb · [Vka (e) − Vka (c)] · Vlb (c) · (B̄ab log B̄kl − B̄kl )

lb

X

πa πb · (Vlb (e) − Vlb (c)) · (B̄ab log B̄al − B̄al )

a

πa πb (Vlb (e) − Vlb (c)) · (B̄ab log B̄al − B̄al ).

abl

Note that for any V (e) ∈ V, it holds that Vbb (e) − 1 = −

P

l6=b Vlb (e).

We can then

continue the calculation from (26):
(27)
g 0 (0; V (e))
i
X
Xh
=2
πa πb · (Vbb (e) − 1) · (B̄ab log B̄ab − B̄ab ) +
πa πb · Vlb (e) · (B̄ab log B̄al − B̄al )
ab

=2

l6=b

XX

πa πb · Vlb (e) · [B̄ab log B̄al − B̄al − B̄ab log B̄ab + B̄ab ]

ab l6=b
(a)

=2

XX

πa πb · Vlb (e) ·

ab l6=b

= − c̃2 ·

XX
b

l6=b

[Vlb (e) ·

XX
−B̄ab (B̄al − B̄ab )2 (b)
≤
−c̃
·
Vlb (e) · (B̄al − B̄ab )2
2
2
2B̃ab
ab l6=b

X
XX
(c)
c̃3
(B̄al − B̄ab )2 ] ≤ −c̃3 ·
Vlb (e) = − kV (e) − V (c)k1 ,
2
a
b

l6=b

where B̃ab is a number between B̄ab and B̄al . To obtain (a), we have used the
second order Taylor expansion of the function B̄ab log x − x around its maxima
x = B̄ab ; (b) is simply due to mina πa > 0, minab B̄ab > 0; (c) holds since B̄ has
no two identical columns. In order to obtain (24), we need to evaluate g 0 (; V (e))

31

COMMUNITY DETECTION WITH NODAL INFORMATION

for small . From (25) and (26), it is straightforward to see
(28)
g 0 (; V (e)) − g 0 (0; V (e))
X
= 2 ·
πa πb · [Vka (e) − Vka (c)] · [Vlb (e) − Vlb (c)] · (B̄ab log rkl − rkl )
klab

|

{z

}

,G1

+

X

πa πb B̄ab ([Vka (e) − Vka (c)]Vlb (c) + [Vlb (e) − Vlb (c)]Vka (c)) · log

|klab

{z

rkl
B̄kl
}

,G2

−

X

πa πb ([Vka (e) − Vka (c)]Vlb (c) + [Vlb (e) − Vlb (c)]Vka (c)) · (rkl − B̄kl ) .

klab

|

{z

}

,G3

We now bound G1 , G2 and G3 in the above equation. For G1 , note that | log rkl | ≤
maxab | log B̄ab | and |rkl | ≤ maxab B̄ab . Therefore, ∃ c̃4 > 0 such that
(29)
|G1 | ≤ c̃4  ·

XX
ka

|Vka (e) − Vka (c)| · |Vlb (e) − Vlb (c)| ≤ 2K 2 c̃4  · kV (e) − V (c)k1 .

bl

Regarding G2 , since Ṽka − Vka (c) = O(), we have
log

rkl
πk πl B̄kl + O()
= O().
= log
B̄kl
πk πl B̄kl + O()

So we can bound G2 :
(30)
|G2 | ≤ O() ·

X
(|Vka (e) − Vka (c)| + |Vlb (e) − Vlb (c)|) = O() · kV (e) − V (c)k1 .
klab

Similar arguments can give us |G3 | = O() · kV (e) − V (c)k1 . This fact combined
with (27), (28), (29) and (30) completes the proof of (24). We now consider any
V (e) such that kV (e) − V (c)k1 ≤ c̃1 . Define V ∗ (e) = V (c) +

V (e)−V (c)
kV (e)−V (c)k1 .

It is

32

H. WENG AND Y. FENG

then straightforward to confirm that V ∗ (e) ∈ V. Hence,
F (T (c), f0 (c)) − F (T (e), f0 (e)) = H(V (c)) − H(V (e))
= g(0; V ∗ (e)) − g(kV (e) − V (c)k1 ; V ∗ (e))
(e)

(d)

= −g 0 (˜
; V ∗ (e)) · kV (e) − V (c)k1 ≥ c̃0 · kV ∗ (e) − V (c)k1 · kV (e) − V (c)k1

= c̃0 · kV (e) − V (c)k1 ,
where (d) is simply by mean value theorem; ˜ is between 0 and kV (e − V (c))k1 ;
(e) holds because of (24). This finishes the proof.

Proof of Theorem 1. It is not hard to check that similar proofs as the ones of
Lemma 1 and Theorem 2 in Bickel et al. (2013) work under NSBM4 . For simplicity, we do not repeat the derivations here. As a result, we can obtain that as
n → ∞,
(31)
√

p

n(β̂ − β0 ) → N (0, I −1 (β0 )),

n2 ρn log

B̂ab
0 = Op (1), 1 ≤ a, b ≤ K.
Bab

Based on (31), we would like to show the strong consistency of ĉ. Define
N (e; {Bab }) =

X h Oab (e)
µn

ab

log

i
Bab na (e)nb (e) − Oab (e)
+
log(1 − Bab ) ,
ρn
µn

K
i
X
1 Xh T
βei xi − log
exp(βkT xi ) .
C(e; β) =
µn
i

k=1

We then easily see that
ĉ = arg max N (e; {B̂ab }) + C(e; β̂).
e∈{1,...,K}n
4

The tail condition on x is used to show that Theorem 1 in Bickel et al. (2013) holds under

NSBM.

33

COMMUNITY DETECTION WITH NODAL INFORMATION

The subsequent proof is aligned with the ideas of proving strong consistency in
Zhao et al. (2012). We first prove that ĉ is weak consistent. Note that
(32)
0
max |N (e; {B̂ab }) − H(T (e), f 0 (e))| ≤ max |N (e; {B̂ab }) − N (e; {Bab
})|
e

e

0
+ max |N (e; {Bab
}) − H(T̂ (e), fˆ(e))|
e

+ max |H(T̂ (e), fˆ(e)) − H(T (e), f 0 (e))|.
e

We aim to bound the three terms on the right hand side of the above inequality. For
the first one, we have
"

max |N (e; {B̂ab }) −
e

0
})|
N (e; {Bab

≤

0 )
B̂ab (1 − Bab
Oab (e)
· log 0
e
µn
Bab (1 − B̂ab )
ab
#
1 − B̂ab
1
.
log
+
0
ρn
1 − Bab

X

According to (31), if we can show maxe

max

Oab (e)
µn

= Op (1), the above inequality will

imply ∃ an → 0 s.t.
0
})| ≤ an ) → 1.
P (max |N (e; {B̂ab }) − N (e; {Bab

(33)

e

For that purpose, we first apply (22) by choosing n = (nρn )−1/3 to have


Okl (e)
P max max
− T̂kl (e) ≤ n → 1.
e
kl
µn
P
Also notice that maxe |Tkl (e)| ≤ ab B̄ab , and
(34)

(35)

max |Tkl (e) − T̂kl (e)| ≤
e

X

P

B̄ab (|π̂a − πa | + |π̂b − πb |) → 0.

ab

Combining (34) and (35) yields the result maxe

Oab (e)
µn

= Op (1). Regarding the

second term on the right hand side of (32), it is straightforward to see that
X
Okl (e)
0
max |N (e; {Bab
}) − H(T̂ (e), fˆ(e))| ≤ max max
− T̂kl (e) ·
| log B̄ab |+
e
e
kl
µn
ab

X log(1 − ρn B̄ab )
Okl (e) X
·
| log(1 − ρn B̄ab )| +
+ B̄ab .
max max
e
kl
µn
ρn
ab

ab

34

H. WENG AND Y. FENG

The fact that maxe maxkl

Okl (e)
µn

= Op (1) and ρn → 0, combined with (34) en-

ables us to conclude from the last inequality: ∃ bn → 0 such that
0
P (max |N (e; {Bab
}) − H(T̂ (e), fˆ(e))| ≤ bn ) → 1.

(36)

e

For the third term on the right hand side of (32), it is easily seen that
max |H(T̂ (e), fˆ(e)) − H(T (e), f 0 (e))| ≤
e

X
ab

X

max |Tab (e) − T̂ab (e)| · | log B̄ab |+
e

ab

max(|fa0 (e) − fˆa (e)| + |fb0 (e) − fˆb (e)|) · B̄ab = O(max |π̂a − πa |).
e

a

Hence there exists cn → 0 such that
(37)

P (max |H(T̂ (e), fˆ(e)) − H(T (e), f 0 (e))| ≤ cn ) → 1.
e

Putting (32), (33), (36) and (37) together, we obtain that ∃ ˜n → 0 so that
(38)

P (max |N (e; {B̂ab }) − H(T (e), f 0 (e))| ≤ ˜n /2) → 1.
e

We now turn to bounding maxe |C(e; β̂)|. We first decompose it as
C(e; β̂) =

i
n h µn
1X T
C(c; β̂) +
(β̂ei − β̂cTi )xi .
µn n
n
i

According to (31) and (59), it is not hard to get

µn
n C(c; β̂)

= Op (1). Also a direct

Cauchy-Schwartz inequality leads to
max
e

1X T
1X
(β̂ei − β̂cTi )xi ≤ 2kβ̂k2 ·
kxi k2 = Op (1).
n
n
i

i

Therefore, we have (choosing ˜n large enough)
P (max |C(e; β̂)| ≤ ˜n /2) → 1.

(39)

e

The two high probability results (38) and (39) together give us
P (max |N (e; {B̂ab }) + C(e; β̂) − H(T (e), f 0 (e))| ≤ ˜n ) → 1.
e

COMMUNITY DETECTION WITH NODAL INFORMATION

35

Based on this result, Lemma 1 implies that ∃ δn = O(˜
n ) s.t
P



max

kV (e)−V (c)k1 >δn


[N (e; {B̂ab }) + C(e; β̂)] < N (c; {B̂ab }) + C(c; β̂) → 1,

which leads to the weak consistency of ĉ. To obtain strong consistency, it suffices
to show
P



max

0<kV (e)−V (c)k1 ≤δn


[N (e; {B̂ab }) + C(e; β̂)] < N (c; {B̂ab }) + C(c; β̂) → 1.

Denote E = {e : 0 < kV (e) − V (c)k1 ≤ δn }. We first bound N (e; {B̂ab }) −
N (c; {B̂ab }).
N (e; {B̂ab }) − N (c; {B̂ab })
0
0
= N (e; {B̂ab }) − N (e; {Bab
}) + N (c; {Bab
}) − N (c; {B̂ab })
{z
}
|
,N1 (e)

+ H(T (e), f0 (e)) − H(T (c), f0 (c))
|
{z
}
,N2 (e)
0
0
}) .
}) − H(T (e), f0 (e)) + H(T (c), f0 (c)) − N (c; {Bab
+ N (e; {Bab
|
{z
}
,N3 (e)

For N1 (e), we have
(40)
max |N1 (e)| ≤ max max
e∈E

e∈E

ab

max max
e∈E

ab

0 )
B̂ab (1 − Bab
|Oab (e) − Oab (c)| X
+
·
log 0
µn
Bab (1 − B̂ab )
ab

|na (e)nb (e) − na (c)nb (c)| X
1 − B̂ab 1
·
log
0 ρ .
2
n
1 − Bab
n

From (31) we see that log

ab

0 )
B̂ab (1−Bab
0 (1−B̂ )
Bab
ab

B̂ab
= op (1), ρ1n · log 1−
= op (1). It is also
1−B 0
ab

straightforward to confirm,
(41) max
ab

2X
|na (e)nb (e) − na (c)nb (c)|
≤
1(ei 6= ci ) ≤ kV (e) − V (c)k1 .
n2
n
i

36

H. WENG AND Y. FENG

Moreover, note that there exists C1 > 0 such that
(42)

max max |T̂kl (e) − T̂kl (c)|
e∈E kl
XX
≤ max max max B̄ab ·
(|Vka (e) − Vka (c)| + |Vlb (e) − Vlb (c)|)
kl

e∈E

ab

ka

lb

≤ C1 · max kV (e) − V (c)k1 .
e∈E

We can then derive from (23) and (42) that

i

h
Oab (e) Oab (c)
(43) P max max
−
− C2 · kV (e) − V (c)k1 ≤ 0 → 1,
ab
e∈E
µn
µn
where C2 > 0 is a constant. Combining (40), (41) and (43) yields
P (max |N1 (e)| ≤ o(1) · max kV (e) − V (c)k1 ) → 1.

(44)

e∈E

e∈E

Regarding N2 (e), Lemma 1 shows that ∃ C3 > 0 s.t.
N2 (e) ≤ −C3 · kV (e) − V (c)k1 ,

(45)

for e ∈ E .

For the last term, we can express N3 (e) as
(46)
N3 (e) =

X  Oab (e)
ab

+

µn

−


Oab (c)
− T̂ab (e) + T̂ab (c) · log B̄ab
µn

X Oab (c) − Oab (e)
ab

µn

· log(1 − B̄ab ρn )

X na (e)nb (e) − na (c)nb (c) 

1
+
· B̄ab +
log(1 − B̄ab ρn )
2
n
ρn
ab
X
+
(T̂ab (e) − Tab (e) − T̂ab (c) + Tab (c)) · log B̄ab
ab

+

X

(fa0 (e)fb0 (e) − fˆa (e)fˆb (e) − fa0 (c)fb0 (c) + fˆa (c)fˆb (c)) · B̄ab .

ab

With a few steps of calculations, it is not hard to obtain
max |T̂ab (e) − Tab (e) − T̂ab (c) + Tab (c)| ≤ kV (e) − V (c)k1 · O(max |π̂a − πa |),
ab

a

COMMUNITY DETECTION WITH NODAL INFORMATION

37

and
max |fa0 (e)fb0 (e) − fˆa (e)fˆb (e) − fa0 (c)fb0 (c) + fˆa (c)fˆb (c)|
ab

≤ kV (e) − V (c)k1 · O(max |π̂a − πa |).
a

We can then use the above results, together with (23), (41) and (43) to bound the
terms on the right hand side of (46). As a result, we are able to show that for any
C4 > 0,
P (max [|N3 (e)| − C4 · kV (e) − V (c)k1 ] ≤ 0) → 1.

(47)

e∈E

The bounds on N1 (e), N2 (e), N3 (e) in (44), (45) and (47) enables us to conclude
that ∃ C5 > 0,
(48) P (max [N (e; {B̂ab }) − N (c; {B̂ab }) + C5 · kV (e) − V (c)k1 ] ≤ 0) → 1.
e∈E

As a next step, we bound C(e; β̂)−C(c; β̂). The result (66) implies that n maxµinkxi k2 =
op (1). This leads to
n K
1 XX T
C(e; β̂) − C(c; β̂) =
β̂k xi (1(ei = k) − 1(ci = k))
µn

≤

i=1 k=1
n
2kβ̂k2 X

µn

kxi k2 · 1(ei 6= ci )

i=1

2kβ̂k2 n maxi kxi k2
·
≤
µn

Pn

i=1

1(ei 6= ci )
n

(a)

= op (1) · kV (e) − V (c)k1 ,

where (a) holds because kβ̂k2 = Op (1), and

1
n

Pn

i=1

1(ei 6= ci ) ≤ 12 kV (e) −

V (c)k1 . Hence we know


C5
(49)
P max [C(e; β̂) − C(c; β̂) −
· kV (e) − V (c)k1 ] ≤ 0 → 1.
e∈E
2
Finally, from the results of (48) and (49), we are able to show that
P (max [N (e; {B̂ab }) + C(e; β̂)] < N (c; {B̂ab }) + C(c; β̂)) → 1.
e∈E

38

H. WENG AND Y. FENG

This completes the proof of strong consistency for ĉ.


Proof of Theorem 2. It is not difficult to verify that the same results of Lemma 3
and Theorem 3 in Bickel et al. (2013) hold under NSBM. We can thus conclude,
√

p
B̌ab
n2 ρn log 0 = Op (1), 1 ≤ a, b ≤ K.
Bab

n(β̌ − β0 ) → N (0, I −1 (β0 )),

The proof for the strong consistency of č follows the same lines of arguments as
for ĉ in Theorem 1. We hence do not repeat the steps here.

Proof of Theorem 3. We first focus on analyzing c̃. Define
N (e) =

i
Oab (e)
1 Xh
Oab (e) log
− Oab (e) ,
µn
na (e)nb (e)ρn
ab

C(e) =

K
X
i
Xh
1
max
βeTi xi − log
exp(βkT xi ) .
µn β∈RKp ,
βK =0

i

k=1

According to (9), c̃ can be equivalently expressed as
c̃ = arg max [N (e) + C(e)].
e∈{1,...,K}n

Under the conditions in part (i), we claim that the following results hold: ∃ C1 , C2 >
0 such that
(A) P (maxe |N (e) + C(e) − F (T (e), f 0 (e))| < C1 (nρn )−1//2 ) → 1, as n →
∞.
(B) for larger n,
min

{e:kV (e)−V

[F (T (c), f 0 (c))−F (T (e), f 0 (e))] ≥ 3C1 (nρn )−1/2 .

3C
(c)k1 ≥ C 1 (nρn )−1/2 }
2

Suppose for now the two results above are correct. Choosing δn =

3C1
−1/2 ,
C2 (nρn )

COMMUNITY DETECTION WITH NODAL INFORMATION

39

we can obtain that
max

[N (e) + C(e) − N (c) − C(c)]

max

[N (e) + C(e) − F (T (e), f 0 (e))]

{e:kV (e)−V (c)k1 ≥δn }

≤

{e:kV (e)−V (c)k1 ≥δn }

+

max

[F (T (e), f 0 (e)) − F (T (c), f 0 (c))]

{e:kV (e)−V (c)k1 ≥δn }

+ [F (T (c), f 0 (c)) − N (c) − C(c)]
≤ 2 max |N (e) + C(e) − F (T (e), f 0 (e))| − 3C1 (nρn )−1/2 ,
e

for large enough n. Hence as n → ∞,


P
max
[N (e) + C(e)] < N (c) + C(c)
{e:kV (e)−V (c)k1 ≥δn }

≥ P (max |N (e) + C(e) − F (T (e), f 0 (e))| < C1 (nρn )−1/2 ) → 1.
e

This would lead to the first result in part (i). Regarding the proof of (A), we follow
similar arguments in the proof of Theorem 4.1 of Zhao et al. (2012). To save space,
we do not detail all the calculations. The key steps are to show that ∃ C3 , C4 , C5 >
0, such that as n → ∞
(50)

P (max max |Oab (e)/µn − T̂ab (e)| < C3 (nρn )−1/2 ) → 1,

(51)

P (max |π̂a − πa | < C4 (nρn )−1/2 ) → 1,

(52)

P (max |C(e)| ≤ C5 (nρn )−1/2 ) → 1.

e

ab

a

e

The result (50) can be obtained from (22) by choosing C3 sufficiently large. The
condition ρn → 0 combined with a direct application of Hoeffding’s inequality
gives (51). And it is straightforward to verify that for any e ∈ {1, 2, . . . , K}n ,
0 ≥ C(e) ≥

K
X
i −n log K
1 Xh T
0ei xi − log
exp(0Tk xi ) =
.
µn
µn
i

k=1

This together with the condition nρn → ∞ yields (52). To prove (B), we first use
Lemma 1:
(53)
F (T (c), f 0 (c)) − F (T (e), f 0 (e)) ≥ C6 · kV (e) − V (c)k1 , if kV (e) − V (c)k1 ≤ C7 ,

40

H. WENG AND Y. FENG

where C6 , C7 are two positive constants. Furthermore, since V (c) is the unique
maximizer of F (T (e), f 0 (e)) as a function of V (e) ∈ V (see the proof of Theorem
3.4 in Zhao et al. (2012)), we have
(54)

[F (T (c), f 0 (c)) − F (T (e), f 0 (e))] > 3C1 (nρn )−1/2 ,

min

{e:kV (e)−V (c)k1 >C7 }

for large n. On the other hand, (53) implies that
(55)
min
{e: C7 ≥kV (e)−V

3C
(c)k1 ≥ C 1 (nρn )−1/2 }
6

[F (T (c), f 0 (c)) − F (T (e), f 0 (e))] ≥ 3C1 (nρn )−1/2 .

Finally (54) and (55) together finishes the proof of (B) with the choice C2 = C6 .
The next step is to prove the strong consistency of c̃ in part (ii). To derive strong
consistency, we also need to show
(56)

P




[N (e) + C(e)] < N (c) + C(c) → 1,

max

{e:0<kV (e)−V (c)k1 <δn }

which requires a refined analysis. Denote E = {e : 0 < kV (e) − V (c)k1 < δn }.
We make use of the existing result: ∃ C7 > 0 such that as n → ∞,

P max[N (e) − N (c) + C7 · kV (e) − V (c)k1 ] < 0 → 1.
e∈E

This is obtained by combining (A.13) in Zhao et al. (2012) and (1.1) in Bickel et al.
(2015). If we can show
(57)


h
i

P max C(e) − C(c) − C7 · kV (e) − V (c)k1 ≤ 0 → 1,
e∈E

the result (56) will be proved by simply putting together the last two high probability arguments. To derive (57), we introduce several notations. Denote
" K
#
K
X

X
T
T
R(β) = E
βk x1(c = k) − log
exp(βk x) ,
k=1

1
Rn (β, e) =
n

n
X
i=1

k=1

"

K
X
k=1

βkT xi 1(ei = k) − log

K
X
k=1

#

exp(βkT xi ) .

COMMUNITY DETECTION WITH NODAL INFORMATION

41

Observe that Rn (β, c) is the sample version of R(β). We further define an Mestimator:
β̂(e) = arg max Rn (β, e).

(58)

β

According to the Convexity Lemma (Lemma 7.75) in Liese and Miescke (2007),
P

since −Rn (β, c) is a convex stochastic process and Rn (β, c) → R(β), we can
obtain
P

|Rn (β, c) − R(β)| → 0,

sup

(59)

kβ−β0 k2 ≤

where  is an arbitrary positive constant. Also note that
sup |Rn (β, e) − R(β)| ≤

sup

(60)

kβ−β0 k2 ≤ e∈E
n

+

sup

|Rn (β, c) − R(β)|

sup
kβ−β0 k2 ≤

sup

kβ−β0 k2 ≤ e∈E

K

1 XX T
βk xi (1(ci = k) − 1(ei = k)) ,
n
i=1 k=1

and
n

K

1 XX T
sup sup
βk xi (1(ci = k) − 1(ei = k))
kβ−β0 k2 ≤ e∈E n

(61)

i=1 k=1

n
1X
≤ C8 · sup
kxi k2 · 1(ci 6= ei )
e∈E n

(a)

i=1

(b)

≤ C8 ·
≤ C9 ·

n
1 X

n
1
n

i=1
n
X

kxi kα2
kxi kα2

1

α

· sup

n
1 X

e∈E

1

α

· (δn )

n
α−1
α

1(ci 6= ei )

 α−1
α

i=1
P

→ 0,

i=1

where C8 , C9 are two positive constants; (a) is by Cauchy-Schwarz inequality and
(b) is due to Hölder’s inequality. Putting (59), (60) and (61) together, we are able
to show
(62)

sup

P

sup |Rn (β, e) − R(β)| → 0.

kβ−β0 k2 ≤ e∈E

Having the uniform convergence of Rn (β, e), we next study the convergence of
β̂(e). However, the uniformity in (62) only holds over a compact set. The result

42

H. WENG AND Y. FENG

may not be applied directly. For this reason, we introduce an ancillary variable:
β̄(e) = αβ̂(e) + (1 − α)β0 , where α =


.
+kβ̂(e)−β0 k2

Clearly,

kβ̄(e) − β0 k2 = αkβ̂(e) − β0 k2 ≤ .
Besides, since Rn (β, e) is a concave function of β and β̂(e) is its maximizer , we
have
Rn (β̄(e), e) ≥ αRn (β̂(e), e) + (1 − α)Rn (β0 , e) ≥ Rn (β0 , e).
Based on the last two inequalities and the fact that β0 is the maximizer of R(β),
we can derive
0 ≤ sup [R(β0 ) − R(β̄(e))]
e∈E

≤ sup [R(β0 ) − Rn (β0 , e)] + sup [Rn (β0 , e) − Rn (β̄(e), e)]
e∈E

e∈E

+ sup [Rn (β̄(e), e) − R(β̄(e))] ≤ 2 sup
e∈E

sup

P

|Rn (β, e) − R(β)| → 0.

e∈E kβ−β0 k2 ≤

Hence,
P

sup [R(β0 ) − R(β̄(e))] → 0.

(63)

e∈E

Moreover, it can be directly verified that the maximizer β0 is unique and isolated.
So ∀ ˜ > 0, ∃ η > 0 such that




P sup kβ̄(e) − β0 k2 > ˜ ≤ P sup [R(β0 ) − R(β̄(e))] > η → 0,
e∈E

e∈E

where the last limit is implies by (63). We thus have shown that supe∈E kβ̄(e) −
P

β0 k2 → 0. It then leads to the consistency of β̂(e):
(c)

sup kβ̂(e) − β0 k2 = sup
e∈E

e∈E

 supe∈E kβ̄(e) − β0 k2 P
kβ̄(e) − β0 k2
≤
→ 0,
 − kβ̄(e) − β0 k2
 − supe∈E kβ̄(e) − β0 k2

COMMUNITY DETECTION WITH NODAL INFORMATION

43

where (c) is due to the definition of β̄(e). We are now in the position to derive (57).
C(e) − C(c)
i
nh
=
Rn (β̂(e), e) − Rn (β̂(c), c)
µn
i
nh
Rn (β̂(e), e) − Rn (β̂(e), c) + Rn (β̂(e), c) − Rn (β̂(c), c)
=
µn
n K
i
nh
1 XX T
≤
Rn (β̂(e), e) − Rn (β̂(e), c) =
β̂k (e)xi · [1(ei = k) − 1(ci = k)]
µn
µn
i=1 k=1

n
n
2 X T
2kβ̂(e)k2 X
≤
|β̂k (e)xi | · 1(ei 6= ci ) ≤
kxi k2 · 1(ei 6= ci )
µn
µn
i=1

i=1

n
nkβ̂(e)k2 · maxi kxi k2
2kβ̂(e)k2 · maxi kxi k2 X
1(ei 6= ci ) ≤
kV (e) − V (c)k1 .
≤
µn
µn
i=1

Therefore, we can obtain
h
i
sup C(e) − C(c) − C7 · kV (e) − V (c)k1

(64)

e∈E

≤ sup
e∈E


 nkβ̂(e)k · max kx k
2
i
i 2
− C7 kV (e) − V (c)k1 .
µn
P

Because supe∈E kβ̂(e) − β0 k2 → 0, it is straightforward to show
P (sup kβ̂(e)k2 ≤ 2kβ0 k2 ) → 1.

(65)

e∈E

According to the condition P (kxk2 > t) = O(e−κ2 t ) and

nρn
log n

→ ∞, we have for

any positive constant c0 ,
nP


n
kxk2 > c0 = O(exp{log n − κ2 c0 nρn }) = o(1),
µn

resulting in
(66)

lim P

n→∞

n

h
n
in
max kxi k2 ≤ c0 = lim 1 − P
kxk2 > c0
n→∞
µn i
µn
= lim e
n→∞

−nP ( µn kxk2 >c0 )
n

= 1.

44

H. WENG AND Y. FENG

Choosing c0 =

C7
2kβ0 k2

and combining (65) and (66) lead to

#

nkβ̂(e)k2 · maxi kxi k2
P sup
− C7 ≤ 0
µn
e∈E
n
C7 
≥ P (sup kβ̂(e)k2 ≤ 2kβ0 k2 ) − P
max kxi k2 >
→ 1.
µn i
2kβ0 k2
e∈E


"

The result (57) can be proved by combining the previous line with (64).
We now turn to deriving the consistency of β̃. We first prove the result in part
(i). We have showed that
P

sup kβ̂(e) − β0 k2 → 0,

P (kV (c̃) − V (c)k1 < δn ) → 1.

e∈E

Hence, for any given  > 0,
P (kβ̃ − β0 k2 < ) = P (kβ̂(c̃) − β0 k2 < )


≥ P sup kβ̂(e) − β0 k2 < , kV (c̃) − V (c)k1 < δn
e∈E




≥ P sup kβ̂(e) − β0 k2 <  − P kV (c̃) − V (c)k1 ≥ δn → 1,
e∈E
P

which shows that kβ̃−β0 k2 → 0. To derive the convergence rate, denote pk (x; β) =
exp(βkT x)
PK
.
T
k=1 exp(βk x)

We may suppress the dependency on x or β when it is clear from

the context. By a first order Taylor expansion of ∇Rn (β̂(c̃), c̃) around β0 , we
have5
(67)
n

n

i=1

i=1

1X
1X
0=
[1(c̃i = k) − pk (xi ; β0 )]xi +
Hk (xi ; β́ k )(β̂(c̃) − β0 ),
n
n

1 ≤ k ≤ K,

where
Hk (x; β) , (pk p1 xxT , . . . , pk pk−1 xxT , −pk (1−pk )xxT , pk pk+1 xxT , . . . , pk pK xxT )
5

For notational simplicity, we have included k = K in the subsequent arguments, though we

should keep in mind that βK ≡ 0.

45

COMMUNITY DETECTION WITH NODAL INFORMATION
P

P

and β́ k is between β̂(c̃) and β0 . Since kβ̂(c̃) − β0 k2 → 06 , kβ́ k − β0 k2 → 0. We
can then apply Theorem 9.4 in Keener (2010) to conclude that
n

1X
P
Hk (xi ; β́ k ) → EHk (x; β0 ),
n

(68)

1 ≤ k ≤ K.

i=1

P
Next we analyze the term n1 ni=1 [1(c̃i = k) − pk (xi ; β0 )]xi in (67). First note
P
that √1n ni=1 [1(ci = k) − pk (xi ; β0 )]xi = Op (1). Also,
n

1X
[1(c̃i = k) − 1(ci = k)]xi
n
≤

1
n

i=1
n
X
i=1

n

2

≤

1X
1(c˜i 6= ci ) · kxi k2
n
i=1

n
 α−1  1 X
1
α
α
·
1(c̃i =
6 ci )
kxi kα2
n
i=1

Therefore, we can obtain
n

(69)

1−α
1X
[1(c̃i = k) − pk (xi ; β0 )]xi = Op ((nρn ) 2α ).
n

i=1

Based on (67), (68) and (69), it is straightforward to derive the second result of part
(i). Regarding the asymptotic normality of β̃ in part (ii), we have already proved
√
that P (c̃ = c) → 1. From standard asymptotic results, we know that n(β̂(c) −
√
d
β0 ) → N (0, I −1 (β0 )). The proof can be finished by showing that n(β̂(c) −
P

β̃) → 0, with one line of arguments as follows.
√
√
P (| n(β̂(c) − β̃)| < ) = P (| n(β̂(c) − β̂(c̃))| < ) ≥ P (c̃ = c) → 1.

Proof of Theorem 4. The proof is motivated by the arguments presented in Guédon
and Vershynin (2016). To make the analysis compact, we will use the concentration
inequality proved there:
(70)
6



n(n − 1)t 
n(n − 1)t2 
P z T [A − E(A | c)]y >
| c ≤ exp −
,
2
16p̄(c) + 8t/3

We recall that β̂(c̃) = β̃, according to our definition (58).

46

H. WENG AND Y. FENG

where p̄(c) =

2
n(n−1)

P

i<j

Var(Aij | c), z, y ∈ {−1, 1}n . The above result is a

direct application of Bernstein’s inequality (See Page 14 in Guédon and Vershynin
(2016)). Let M (c) ∈ Rn×K with Mik (c) = 1 if ci = k; Mik (c) = 0 otherwise;
Z̄(c) = M (c)M T (c) and
K
n
n
X
X
X
 o
n×n
Mc = Z ∈ R
: Z  0, 0 ≤ Zij ≤ 1,
Zij =
1(ci = k) 2 .
ij

k=1

i=1

It is straightforward to verify that E(A | c) = ρn M (c)B̄M T (c) and Z̄ij (c) = 1 if
and only if ci = cj .
Before proceeding, define Sn (Z) = hA + γn XX T , Zi and S(Z) = hE(A |
c) + γn E(X | c)E(X T | c), Zi. Then, by the definition in (10), we have Ẑ =
arg maxZ∈Mc Sn (Z). Hence for any Z ∈ Mc , we have
(71)
S(Z̄(c)) − S(Z)
= hE(A | c) + γn E(X | c)E(X T | c), Z̄(c) − Zi
X
=
(ρn B̄ci cj + γn E(xTi | ci )E(xj | cj )) · (Z̄ij (c) − Zij ) · 1(ci = cj )
ij

−

X

(ρn B̄ci cj + γn E(xTi | ci )E(xj | cj )) · (Zij − Z̄ij (c)) · 1(ci 6= cj )

ij
(a)

≥U·

X

(Z̄ij (c) − Zij ) · 1(ci = cj ) − L ·

ij
(b)

≥

X
(Zij − Z̄ij (c)) · 1(ci 6= cj )
ij

U −L
· kZ̄(c) − Zk1 ,
2

where U = min1≤k≤K {ρn B̄kk +γn E(xT | c = k)E(x | c = k)}, L = maxa6=b {ρn B̄ab +
γn E(xT | c = a)E(x | c = b)}; in (a) we have used the fact that Z̄ij (c) ≥ Zij if
P
ci = cj and Z̄ij (c) ≤ Zij otherwise; (b) holds because ij (Z̄ij (c) − Zij ) · 1(ci =
P
cj ) = ij (Zij − Z̄ij (c)) · 1(ci 6= cj ) = 12 kZ̄(c) − Zk1 . Note that the conditions
we assumed imply U > L > 0 when n is large enough. In addition, we observe

COMMUNITY DETECTION WITH NODAL INFORMATION

47

that
(72)
S(Z̄(c)) − S(Ẑ) ≤ [S(Z̄(c)) − Sn (Z̄(c))] + [Sn (Z̄(c)) − Sn (Ẑ)] + [Sn (Ẑ) − S(Ẑ)]
≤ 2 max |Sn (Z) − S(Z)|,
Z∈Mc

since Sn (Z̄(c)) − Sn (Ẑ) ≤ 0.
Moreover, by applying the Grothendieck’s inequality described in Section 3.1
of (Guédon and Vershynin, 2016), we can have
max |Sn (Z) − S(Z)|

(73)

Z∈Mc

= max |hA − E(A | c) + γn XX T − γn E(X | c)E(X T | c), Zi|
Z∈Mc

≤2

max

z,y∈{−1,1}n

z T [A − E(A | c) + γn XX T − γn E(X | c)E(X T | c)]y.

Combining (71), (72), and (73), it is not hard to obtain
kZ̄(c) − Ẑk1 ≤

(74)

8
·
max z T [A − E(A | c) + γn XX T
U − L z,y∈{−1,1}n
− γn E(X | c)E(X T | c)]y.

We now bound the term on the right hand side of the above inequality. Note that
for sufficiently large n
p̄(c) =

X
X
2
minab B̄ab
1
min B̄ab
ρn =
· ρn .
ρn B̄ci cj (1 − ρn B̄ci cj ) ≥
n(n − 1)
n(n − 1) ab
2
i<j

i<j

2/3

Therefore, according to (70), choosing t = n−1/3 ρn yields that there exists C1 >
0 such that when n is large



P z T [A − E(A | c)]y > n5/3 ρn2/3 | c ≤ exp − C1 n · (nρn )1/3 ,
for any z, y ∈ {−1, 1}n . Using the union bound, we can conclude
(75)

P

max

z,y∈{−1,1}n



z T [A − E(A | c)]y > n5/3 ρ2/3
≤ 22n · exp − C1 n · (nρn )1/3 → 0.
n

48

H. WENG AND Y. FENG

Regarding the bound on γn XX T − γn E(X | c)E(X T | c), we first denote the i-th
column of X by di . Since di | c has independent sub-Gaussian elements7 , we can
apply Hoeffding’s inequality to obtain

2
P |y T (di − E(di | c))| > nt | c ≤ 2e−C2 nt ,
where C2 > 0 is a constant and y ∈ {1, −1}n . Note that we can choose C2 small
enough to guarantee that the above inequality holds for all c. Accordingly, we have

2
(76)
P |y T (di − E(di | c))| > nt ≤ 2e−C2 nt .
This implies the following bound for large t, y, z ∈ {−1, 1}n ,
P (|z T [di dTi − E(di | c)E(dTi | c)]y| > n2 t2 )

(77)

≤ P (|z T (di − E(di | c)) · y T (di − E(di | c))| > n2 t2 /3)
+ P (|z T (di − E(di | c)) · y T E(di | c)| > n2 t2 /3)
+ P (|z T E(di | c) · y T (di − E(di | c))| > n2 t2 /3)
(c)

2

4

4

2

≤ 4e−C3 nt + 2e−C4 nt + 2e−C4 nt ≤ 8e−C5 nt .

where C3 , C4 , C5 are positive constants; we have used (76) and the fact that
maxy∈{−1,1}n |y T E(di | c)| = O(n) to obtain (c). Applying the inequality (77)
enables us to conclude that
P (γn |z T (XX T − E(X | c)E(X T | c))y| > n2 t2 )
2 /(pγ

≤ p · P (|z T [di dTi − E(di | c)E(dTi | c)]y| > n2 t2 /(pγn )) ≤ 8pe−C5 nt

n)

.

√
Choosing t = C6 γn and using the union bound, we then have


(78)
P
max γn |z T (XX T − E(X | c)E(X T | c))y| > C6 n2 γn
z,y∈{−1,1}n

2

≤ 22n · 8pe−C5 C6 np

−1

→ 0,

where the last limit holds for sufficiently large C6 . Putting (74), (75) and (78)
together gives us that
2/3

kZ̄(c) − Ẑk1
8n−1/3 ρn + 8C6 γn
≤
,
n2
U −L
7

This can be directly shown from the condition that kxk2 is sub-Gaussian.

49

COMMUNITY DETECTION WITH NODAL INFORMATION

holds with probability approaching 1. Since U − L is of order ρn and γn =
o(ρn ), nρn → ∞, we obtain that

kZ̄(c)−Ẑk1 P
→
n2

0.

Finally, we use similar arguments as in Rohe et al. (2011) to analyze the Kmeans step and show the mis-classification rate vanishes. Denote the K centroids
output from K-means run on Ẑ by µ̂1 , . . . , µ̂K and on Z̄(c) by µ1 , . . . , µK ; Ĉ =
(µ̂c̄1 , . . . , µ̂c̄n )T ; C = {C ∈ Rn×n : C has exactly K non-identical rows}. We
then know Z̄(c) = (µc1 , . . . , µcn )T and also
Ĉ = arg min kC − Ẑk2F .

(79)

C∈C

Note that the elements of µi are either 0 or 1 and mutually orthogonal. Since the
size of each community is proportional to n, we have with high probability, there
√
exists a constant C6 > 0 such that kµi − µj k2 ≥ C6 n for 1 ≤ i 6= j ≤ K.
We now define a set of nodes that are “incorrectly” identified by the K-means:
√
N = {i ∈ {1, 2, . . . , n} : kµ̂c̄i − µci k2 ≥ C6 n/2}.
Hence, we have
|N |
4 X
1 X
4
≤
kµ̂c̄i − µci k22
kµ̂c̄i − µci k22 · 2 ≤ 2 2
n
n
C6 n
C6 n i
i∈N

=

(d) 16
4
16
2
k
Ĉ
−
Z̄(c)k
≤ 2 2 kẐ − Z̄(c)k2F ≤ 2 2 kẐ − Z̄(c)k1 ,
F
2
2
C6 n
C6 n
C6 n

where (d) is due to (79). This combined with the fact that
that

|N | P
n →

kZ̄(c)−Ẑk1 P
→
n2

0 shows

0. On the other hand, for any i ∈
/ N , it is straightforward to verify that

for any 1 ≤ k 6= ci ≤ K, we have
kµ̂c̄i − µk k2 ≥ kµci − µk k2 − kµ̂c̄i − µci k2
√
√
√
≥ C6 n − C6 n/2 = C6 n/2 > kµ̂c̄i − µci k2 .
We construct a new community assignment:
c̄¯i = arg min kµ̂c̄i − µa k2 .
a∈{1,...,K}

50

H. WENG AND Y. FENG

Note that the communities that {c̄¯i } represents might be different from {c̄i }, since
{c̄¯i } could cover less than K communities. We can view {c̄i } as a possibly finer
partition of {c̄¯i }. Clearly, for any i ∈
/ N , c̄¯i = ci . We thus have
(80)

1X
1
P
1(c̄¯i 6= ci ) ≤ |N | →
0.
n
n
i

This implies that {c̄¯i } = {c̄i } with probability approaching 1. Otherwise, the estimate {c̄¯i } yields less than K communities with non-vanishing probability, which
is impossible for {c̄¯i } to achieve the result in (80). This completes the proof.

REFERENCES
A BBE , E. and S ANDON , C. (2015). Community detection in general stochastic block models: fundamental limits and efficient recovery algorithms. arXiv:1503.00609.
A IROLDI , E. M., B LEI , D. M., F IENBERG , S. E. and X ING , E. P. (2009). Mixed membership
stochastic blockmodels. In Advances in Neural Information Processing Systems. 33–40.
A KOGLU , L., T ONG , H., M EEDER , B. and FALOUTSOS , C. (2012). Pics: Parameter-free identification of cohesive subgroups in large attributed graphs. In SDM. Citeseer, 439–450.
A MINI , A. A., C HEN , A., B ICKEL , P. J. and L EVINA , E. (2013). Pseudo-likelihood methods for
community detection in large sparse networks. The Annals of Statistics, 41 2097–2122.
A MINI , A. A. and L EVINA , E. (2014).

On semidefinite relaxations for the block model.

arXiv:1406.5647.
A NA , L. and JAIN , A. K. (2003). Robust data clustering. In Computer Vision and Pattern Recognition, 2003. Proceedings. 2003 IEEE Computer Society Conference on, vol. 2. IEEE, II–128.
A NANDKUMAR , A., G E , R., H SU , D. and K AKADE , S. M. (2014). A tensor approach to learning
mixed membership community models. The Journal of Machine Learning Research, 15 2239–
2312.
B ICKEL , P., C HOI , D., C HANG , X. and Z HANG , H. (2013). Asymptotic normality of maximum
likelihood and its variational approximation for stochastic blockmodels. The Annals of Statistics,
41 1922–1943.
B ICKEL , P. J. and C HEN , A. (2009). A nonparametric view of network models and newman–girvan
and other modularities. Proceedings of the National Academy of Sciences, 106 21068–21073.
B ICKEL , P. J., C HEN , A., Z HAO , Y., L EVINA , E. and Z HU , J. (2015). Correction to the proof of
consistency of community detection. The Annals of Statistics, 43 462–466.
B INKIEWICZ , N., VOGELSTEIN , J. T. and ROHE , K. (2014). Covariate assisted spectral clustering.
arXiv preprint arXiv:1411.2158.

COMMUNITY DETECTION WITH NODAL INFORMATION

51

B OYD , S., PARIKH , N., C HU , E., P ELEATO , B. and E CKSTEIN , J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and
Trends R in Machine Learning, 3 1–122.
C AI , T. T. and L I , X. (2015). Robust and computationally feasible community detection in the
presence of arbitrary outlier nodes. The Annals of Statistics, 43 1027–1059.
C ELISSE , A., DAUDIN , J.-J. and P IERRE , L. (2012). Consistency of maximum-likelihood and
variational estimators in the stochastic block model. Electronic Journal of Statistics, 6 1847–
1899.
C HANG , J. and B LEI , D. M. (2010). Hierarchical relational models for document networks. The
Annals of Applied Statistics 124–150.
C HEN , Y., L I , X. and X U , J. (2015). Convexified modularity maximization for degree-corrected
stochastic block models. arXiv:1512.08425.
C HEN , Y., S ANGHAVI , S. and X U , H. (2012). Clustering sparse graphs. In Advances in neural
information processing systems. 2204–2212.
C HOI , D. S., W OLFE , P. J. and A IROLDI , E. M. (2012). Stochastic blockmodels with a growing
number of classes. Biometrika asr053.
C ROSS , R. L. and PARKER , A. (2004). The hidden power of social networks: Understanding how
work really gets done in organizations. Harvard Business Review Press.
DASGUPTA , A., H OPCROFT, J. E. and M C S HERRY, F. (2004). Spectral analysis of random graphs
with skewed degree distributions. In Foundations of Computer Science, 2004. Proceedings. 45th
Annual IEEE Symposium on. IEEE, 602–610.
DAUDIN , J.-J., P ICARD , F. and ROBIN , S. (2008). A mixture model for random graphs. Statistics
and computing, 18 173–183.
D ECELLE , A., K RZAKALA , F., M OORE , C. and Z DEBOROV Á , L. (2011). Asymptotic analysis
of the stochastic block model for modular networks and its algorithmic applications. Physical
Review E, 84 066106.
D EMPSTER , A. P., L AIRD , N. M. and RUBIN , D. B. (1977). Maximum likelihood from incomplete
data via the em algorithm. Journal of the royal statistical society. Series B (methodological) 1–38.
F ORTUNATO , S. (2010). Community detection in graphs. Physics Reports, 486 75–174.
G AO , C., M A , Z., Z HANG , A. Y. and Z HOU , H. H. (2015). Achieving optimal misclassification
proportion in stochastic block model. arXiv:1505.03772.
G U ÉDON , O. and V ERSHYNIN , R. (2016).

Community detection in sparse networks via

grothendiecks inequality. Probability Theory and Related Fields, 165 1025–1049.
H OLLAND , P. W., L ASKEY, K. B. and L EINHARDT, S. (1983). Stochastic blockmodels: First steps.
Social networks, 5 109–137.
H UBERT, L. and A RABIE , P. (1985). Comparing partitions. Journal of classification, 2 193–218.
J IN , J. (2015). Fast community detection by score. The Annals of Statistics, 43 57–89.
J ORDAN , M. I., G HAHRAMANI , Z., JAAKKOLA , T. S. and S AUL , L. K. (1999). An introduction to

52

H. WENG AND Y. FENG
variational methods for graphical models. Machine learning, 37 183–233.

J OSEPH , A. and Y U , B. (2013). Impact of regularization on spectral clustering. arXiv:1312.1733.
K ARRER , B. and N EWMAN , M. E. (2011). Stochastic blockmodels and community structure in
networks. Physical Review E, 83 016107.
K EENER , R. W. (2010). Theoretical Statistics: Topics for a Core Course (Springer Texts in Statistics), vol. 1. Springer.
K RZAKALA , F., M OORE , C., M OSSEL , E., N EEMAN , J., S LY, A., Z DEBOROV Á , L. and Z HANG , P.
(2013). Spectral redemption in clustering sparse networks. Proceedings of the National Academy
of Sciences, 110 20935–20940.
L E , C. M. and L EVINA , E. (2015). Estimating the number of communities in networks by spectral
methods. arXiv preprint arXiv:1507.00827.
L EI , J. (2016). A goodness-of-fit test for stochastic block models. The Annals of Statistics, 44
401–424.
L EI , J. and R INALDO , A. (2014). Consistency of spectral clustering in stochastic block models. The
Annals of Statistics, 43 215–237.
L IESE , F. and M IESCKE , K.-J. (2007). Statistical decision theory. In Statistical Decision Theory.
Springer, 1–52.
M ONTANARI , A. and S EN , S. (2015). Semidefinite programs on sparse random graphs and their
application to community detection. arXiv:1504.05910.
NALLAPATI , R. and C OHEN , W. W. (2008). Link-plsa-lda: A new unsupervised model for topics
and influence of blogs. In ICWSM.
N EWMAN , M. and C LAUSET, A. (2015).

Structure and inference in annotated networks.

arXiv:1507.04001.
N EWMAN , M. E. (2003). The structure and function of complex networks. SIAM review, 45 167–
256.
N EWMAN , M. E. (2006). Modularity and community structure in networks. Proceedings of the
National Academy of Sciences, 103 8577–8582.
Q IN , T. and ROHE , K. (2013). Regularized spectral clustering under the degree-corrected stochastic
blockmodel. In Advances in Neural Information Processing Systems. 3120–3128.
ROHE , K., C HATTERJEE , S. and Y U , B. (2011). Spectral clustering and the high-dimensional
stochastic blockmodel. The Annals of Statistics, 39 1878–1915.
RUAN , Y., F UHRY, D. and PARTHASARATHY, S. (2013). Efficient community detection in large
networks using content and links. In Proceedings of the 22nd international conference on world
wide web. International World Wide Web Conferences Steering Committee, 1089–1098.
S AADE , A., K RZAKALA , F. and Z DEBOROV Á , L. (2014). Spectral clustering of graphs with the
bethe hessian. In Advances in Neural Information Processing Systems. 406–414.
S ALDANA , D. F., Y U , Y. and F ENG , Y. (2015). How many communities are there? Journal of
Computational and Graphical Statistics. To appear.

COMMUNITY DETECTION WITH NODAL INFORMATION

53

S TEINHAEUSER , K. and C HAWLA , N. V. (2010). Identifying and evaluating community structure
in complex networks. Pattern Recognition Letters, 31 413–421.
T ÜT ÜNC Ü , R. H., T OH , K. C. and T ODD , M. J. (2003). Solving semidefinite-quadratic-linear
programs using sdpt3. Mathematical programming, 95 189–217.
WANG , Y. and B ICKEL , P. J. (2015). Likelihood-based model selection for stochastic block models.
arXiv preprint arXiv:1502.02069.
YANG , J., M C AULEY, J. and L ESKOVEC , J. (2013). Community detection in networks with node
attributes. In Data Mining (ICDM), 2013 IEEE 13th International Conference on. IEEE, 1151–
1156.
Z HANG , A. Y. and Z HOU , H. H. (2015). Minimax rates of community detection in stochastic block
models. arXiv preprint arXiv:1507.05313.
Z HANG , Y., L EVINA , E. and Z HU , J. (2013). Community detection in networks with node features.
In Advances in Neural Information Processing Systems.
Z HANG , Y., L EVINA , E. and Z HU , J. (2014). Detecting overlapping communities in networks with
spectral methods. arXiv:1412.3432.
Z HAO , Y., L EVINA , E. and Z HU , J. (2012). Consistency of community detection in networks under
degree-corrected stochastic block models. The Annals of Statistics, 40 2266–2292.
D EPARTMENT OF S TATISTICS ,
C OLUMBIA U NIVERSITY,
N EW YORK , NY, 10027.
E- MAIL : hw2375@columbia.edu
yang.feng@columbia.edu

