arXiv:1709.00407v2 [stat.ML] 23 Nov 2017

Estimating Mixed Memberships with Sharp Eigenvector
Deviations
Xueyu Mao∗, Purnamrita Sarkar†, Deepayan Chakrabarti‡
The University of Texas at Austin

Abstract
We consider the problem of estimating overlapping community memberships.
Existing provable algorithms for this problem either make strong assumptions
about the population [33, 16], or are too computationally expensive [3, 15]. We
work under the popular Mixed Membership Stochastic Blockmodel (MMSB) [2].
Using the inherent geometry of this model, we link the inference of overlapping
communities to the problem of ﬁnding corners in a noisy rotated and scaled
simplex, for which consistent algorithms exist [12]. We use this as a building block for our algorithm to infer the community memberships of each node.
Furthermore, we prove that each node’s soft membership vector converges to
its population counterpart. To our knowledge, this is the ﬁrst work to obtain
rate of convergence for community membership vectors of each node, in contrast to previous work which obtain convergence results for memberships of all
nodes as a whole. As a byproduct of our analysis, we derive sharp row-wise
eigenvector deviation bounds, and provide a cleaning step that improves the
performance signiﬁcantly for sparse networks. We also propose both necessary
and suﬃcient conditions for identiﬁability of the model, while existing methods typically present suﬃcient conditions. The empirical performance of our
method is shown using simulated and real datasets scaling up to 100,000 nodes.

1

Introduction

In most real-world networks, a node belongs to multiple communities. In an university, professors have joint appointments to multiple departments; a movie like
“Dirty Harry” in the Netﬂix recommendation network belongs to action, thriller,
and the drama genre according to Google; in a book recommendation network like
goodreads.com, “To Kill a Mockingbird” can be classiﬁed as a classic, historical
ﬁction, young-adult ﬁction, etc. The goal of community detection is to consistently
infer each node’s community memberships from just the network structure.
∗

Department of Computer Science. Email: xmao@cs.utexas.edu
Department of Statistics and Data Sciences. Email: purna.sarkar@austin.utexas.edu
‡
Department of Information, Risk, and Operations Management. Email: deepay@utexas.edu

†

1

A well-studied variant of this problem assumes that each node belongs to a single
community. For instance, under the Stochastic Blockmodel [14], the probability of
a link between two nodes depends only on their respective communities. Thus,
provably consistent inference under the Stochastic Blockmodel involves ﬁnding the
unknown cluster membership of each node (see [19]) and these are not immediately
applicable for the general problem where a node may belong to multiple communities
to diﬀerent degrees.
In this paper, we work with the popular Mixed Membership Stochastic Blockmodel (MMSB) [2]. This generalizes the Stochastic Blockmodel by letting each node
i have diﬀerent degrees of membership in all communities. In particular, each node
i is associated with a community membership vector θi ∈ RK (θi ≥ 0, kθi k1 = 1),
drawn from a Dirichlet prior. The model for generating the symmetric adjacency
matrix is as follows:
θi ∼ Dirichlet(α)
T

P := ρΘBΘ

α ∈ RK
+,

i ∈ [n]

Aij = Aji ∼ Bernoulli(Pij ) i, j ∈ [n]

For identiﬁability we assume maxij Bij = 1. For assortative community structures,
B has higher values on its diagonal as compared to the oﬀ-diagonal, implying that
two nodes with a high membership in the same community are more likely to form an
edge. In dissortative settings, (bipartite graphs being an extreme case) oﬀ-diagonal
elements are larger than diagonal elements. We allow both ρ and the smallest
singular value of B, denoted by λ∗ (B), to go to zero. Here ρ controls the overall
edge density of the model, whereas λ∗ (B) can be thought of as a stand-in for the
separation between blocks.
P
The quantity α0 = K
a=1 αa controls the level of overlap between members of
diﬀerent communities: as α0 → 0, the MMSB model degenerates to the Stochastic
Blockmodel [14]. The goal of community detection under the MMSB model is to
recover Θ (matrix with the ith row as θiT ) and B from the observed adjacency matrix
A.
Prior work on this problem include MCMC [2] and computationally eﬃcient
variational approximation methods [13] (SVI) which do not have any guarantees of
consistency. A notable family of algorithms that has been shown to be theoretically
consistent uses tensor-based methods [3, 15]. However, these are typically hard to
implement, and has a computational complexity is O(n2 K), which can be prohibitive
for large graphs.
Recently Mao et al. [20] have proposed a provably consistent geometric algorithm
(GeoNMF) for MMSB with diagonal B and α = α0 1K /K. However the guarantees
only work in the dense regime where average degree grows faster than the square
root of n. Other interesting network models for overlapping communities and nonnegative matrix factorization style inference methods which do not have theoretical
guarantees include [7, 30, 31, 25].

2

Zhang et al. [33] propose a provably consistent spectral algorithm (OCCAM) for
a related but diﬀerent model. The authors assume that each community has some
“pure” nodes (which only belong to that community). This is a common assumption
for non-negative matrix factorization methods [4, 20]. The authors also assume
that B is positive semideﬁnite and full rank with equal diagonal entries. Other
assumptions ensure that the k-medians loss function on θi attains its minimum at
the locations of the pure nodes and there is a curvature around this minimum. This
condition is typically hard to check.
We recently became aware of concurrent work [16] which uses a convex hull
ﬁnding algorithm for a degree corrected variant of the MMSB model. The authors
assume that every cluster has Θ(n) pure nodes. Furthermore, they assume that
either the non-pure nodes are relatively few, or when there are many such nodes,
the θi ’s for these nodes can be grouped into a constant number of well separated
partitions. In contrast, we require that each cluster has at least one pure node.
Kaufmann et al. [18] propose a combinatorial algorithm (SAAC) for detecting
overlapping communities for a related model. Consistency is proven for the global
optimum of their loss function, ﬁnding which is computationally intractable. The
authors propose a nonconvex alternating minimization method instead; however,
no consistency analysis is provided for this more scalable and hence more viable
alternative. Other notable examples of related but diﬀerent models include [26],
who show consistency when the overlap between clusters is small.
In light of these existing works, we present an algorithm to infer the parameters
for a full-rank B (which need not be positive deﬁnite) and prove its theoretical
consistency. Our contributions are as follows.
Identifiability: We present both necessary and suﬃcient conditions for identiﬁability of the MMSB model in Section 2. Speciﬁcally, if each of the K communities
has at least one “pure” node, then a full-rank B is suﬃcient for identiﬁability. Surprisingly, even a B of rank K − 1 can be suﬃcient, as long as no row of B is an
aﬃne combination of the other rows of B. On the other hand, if the entries of ρB
are strictly between 0 and 1, then pure nodes are necessary for identiﬁability. Our
suﬃcient condition is more general than that in [33]. To our knowledge, we are the
ﬁrst to report both necessary and suﬃcient conditions for identiﬁability under the
MMSB model.
Recovery algorithm: It is not hard to see that for the MMSB model, the population eigenvectors (i.e., eigenvectors of the matrix P) form a rotated and scaled
simplex. We present an algorithm called SPACL, which re-purposes an existing
algorithm [12] for detecting corners in a rotated and scaled simplex to ﬁnd pure
nodes, and then uses these to infer the parameters Θ and B. It also includes a
novel preprocessing step that improves performance in sparse settings. The main
compute-intensive parts of the algorithm are a) top-K eigen-decomposition of A,
b) calculating k-nearest neighbors of a point for preprocessing. There are highly
optimized algorithms and data structures for both of these steps [8, 24, 9].
3

Strong consistency: All of the existing works on MMSB type models show consistency in terms of the relative deviation or correlation of Θ̂ as a whole w.r.t. the
truth [33, 18, 15, 16], or columns of Θ̂ [3] (soft memberships of all nodes to a particular community). This can be thought of as a variant of weak consistency in
clustering, where one proves that the fraction of misclassiﬁed nodes go to zero as
the number of nodes grow. In contrast, for any node, we obtain rates of convergence
of the cluster membership vector θ̂i to its population counterpart. This is a variant
of strong consistency in clustering, where one shows that every node is correctly
classiﬁed in the limit. To our knowledge this is the ﬁrst work to establish strong
consistency of an estimation algorithm for MMSB models.
The main diﬃculty is that we need row-wise error bounds of empirical eigenvectors, while existing results from random matrix theory (Davis Kahan Theorem [32],
etc.) are typically on the deviation of the eigenvectors (or eigen-spaces) as a whole.
A straightforward application of these bounds give suboptimal convergence rates on
the parameters, which work when the average degree grows faster than square root
of n. To avoid this diﬃculty, we derive a novel result on row-wise error bounds for
the deviation of empirical eigenvectors from their suitably transformed population
counterparts.
Proof technique: First, we use Kato’s technique of equating eigenspaces with
certain contour integrals of the resolvent of a matrix [17]. Second, we show that the
population eigenvectors are delocalizedpwith high probability, that is, no one element
of any eigenvector has more than Θ( ν/n) mass, where ν := α0 / mink αk . Note
that this is a high probability statement, unlike in Stochastic Blockmodels, because
the population eigenvectors and the θi are random quantities that depend on draws
from a Dirichlet distribution. The third ingredient is a path counting argument
from recent random matrix theory literature [11] which can be easily adapted to
our setting. This analysis is fairly general and can be applied to other random
graph models like Stochastic Blockmodels. These row-wise eigenvector bounds also
enable us to obtain convergence results akin to strong consistency in clustering as
discussed below. The main theorems are presented in Section 3 with details deferred
to Section 5.
Empirical validation: In Section 4, we compare SPACL with OCCAM, variational methods, SAAC and existing non-negative matrix factorization algorithms
(GeoNMF, BSNMF) on both simulated and large real world author-author and
paper-author networks with up-to 100,000 nodes. These comparisons show that
SPACL is better able to recover the model parameters, while being able to process
networks with O(105 ) nodes in tens of seconds. We also present both simulated and
real data examples where SPACL succeeds whereas other methods like symmetric
non-negative matrix factorization and OCCAM are inapplicable since they assume
that B is positive deﬁnite.

4

Table 1: Table of notations. K leading eigenvectors of a matrix correspond to K
largest eigenvalues in magnitude.
n

Number of nodes

ρB ∈

[0, 1]K×K

K

Community link probabilities (B =

BT )

α∈

Number of communities
RK×1
+

Dirichlet prior parameters

n×K
Θ ∈ R+

Fractional community memberships

α0

αmin (αmax )

mini∈[K] αi (maxi∈[K] αi )

ν

P

A

Adjacency matrix

P

ρΘBΘT

ρ

Upper bound on Pij

Im

i αi

α0 /αmin

Rn×K

m × m identity matrix

Ê

Diagonal matrix of eigenvalues of A

E

Diagonal matrix of eigenvalues of P

V ∈ Rn×K

True K pure node index rows of V

λK (M)

K th largest eigenvalue of M

Estimated K pure node index rows of V

λ∗ (M)

K th largest singular value of M

Ô ∈ RK×K

Orthogonal matrix

λi

ith largest eigenvalue of P

Π∈

Permutation matrix

λ̂i

ith largest eigenvalue of A

All ones vector of length m

ei

ei (j) = 1(i == j)

VP ∈

RK×K

Vp ∈ RK×K

1m

2

{0, 1}K×K

V̂ ∈

K leading eigenvectors of A
K leading vigenvectors of P

Notations, Identifiability and Algorithms

Before presenting our results on identiﬁability we introduce some notations and
assumptions.
q

n

27 log n , nρ) λ∗ (B)
α0
≤
,
≥
Assumption 2.1. We assume that, ν :=
αmin
2(1 + α0 )
ν
8(1 + α0 )(log n)ξ
for some constant ξ > 1, and α0 = O(1).
√
nρ

min(

The condition λ∗ (B)/ν immediately implies ρn = Ω((log n)2ξ ), since λ∗ (B) ≤
kBk ≤ K ≤ ν, and α0 is a constant. This condition speciﬁes that our algorithm
works in the semi-dense degree regime, where the average degree grows faster than
polylog of n. This is the most common regime where most consistency results on
network clustering are shown [3, 19, 21] etc. If the network is dense, α0 = O(1),
and λ∗ (B) = Θ(1), since ν ≥ K, the ﬁrst and second conditions essentially states
√
that K = o( n), allowing K to grow with n. Finally, note that in the simple case
of ρB = (pn − qn )IK + qn 1K 1TK , we have ρλ∗ (B) = (pn − qn ). Thus the second
condition can be interpreted as a lower bound
 on cluster separation, since in the
p√
ν
n −qn
simple case mentioned above pn = Ω̃ √n .
For any matrix M, we will use Mi , M1:i,: , and M1:i,1:j to denote the ith row, the
submatrix formed by the ﬁrst i rows, and the submatrix formed by the ﬁrst i rows
and ﬁrst j columns respectively. We will use kMk and kMkF to respectively denote
the operator and Frobenius norms of a n × n matrix, kvk to denote the Euclidean
5

norm of a vector v. We use Õ and Ω̃ to denote upper and lower bounds up to
logarithmic factors. Finally we present a consolidated list of notations in Table 1.
We shall now provide necessary and suﬃcient conditions for the identiﬁability of
the MMSB model, and present our algorithm for inferring its parameters from the
observed network.

2.1

Identifiability

Prior work showed identiﬁability for a related model if B is positive-deﬁnite and has
all ones on the diagonal [33]. We will establish suﬃciency conditions under a weaker
condition. In addition, we will also show necessary conditions for identiﬁability. To
our knowledge this is the ﬁrst work that constructs necessary as well as suﬃcient
conditions for the identiﬁability of MMSB.
Deﬁne a pure node as a node which belongs to exactly one community. All nodes
in a Stochastic Blockmodel are pure nodes, since every node belongs to exactly one
community. Deﬁne a “mixed” node as a node m such that θmj > 0 for all j ∈ [K].
Theorem 2.1. Suppose there are K communities, with at least one pure node for
each community. Then,
(a) If rank(B) = K, then the MMSB model is identifiable up to a permutation.
(b) If rank(B) = K − 1, and no row of B is an affine combination of the other
rows of B, then the MMSB model is identifiable up to a permutation.
(c) In any other case, if there exists a mixed node, then the model is not identifiable.
Our proof links the MMSB parameters Θ and B to the eigen-decomposition of
the probability matrix P, and then exploits its geometric structure. Speciﬁcally, we
will show that the eigenvector row corresponding to any node lies inside a polytope
whose vertices correspond to pure nodes. When B is full rank, the polytope has K
linearly independent vertices, and the community memberships θi of each node i are
ﬁxed by the position of its eigenvector row with respect to these vertices. This proves
part (a). When B is rank-deﬁcient, the points corresponding to the pure nodes are
linearly dependent. However, under the conditions of part (b), the constraints on
Θ are shown to make the model identiﬁable. In other cases, we construct a new Θ′
that still yields the same probability matrix P. This proves part (c).
The next theorem shows that the existence of pure nodes is necessary in most
practical scenarios.
Theorem 2.2. Suppose that ρBij ∈ (0, 1) for all i, j ∈ [K]. MMSB is identifiable up
to a permutation only if there is at least one pure node for each of the K communities.

6

2.2

Algorithms

We will consider the problem of inferring the MMSB parameters under the following
condition.
Assumption 2.2. B ∈ RK×K is full rank, and there is at least one pure node for
each of the K communities.
This is suﬃcient for identiﬁability (Theorem 2.1). Note that this is weaker than
assuming a positive-deﬁnite B, as in [33, 20] or Θ(n) pure nodes, as in [16]. The
addition of one pure node per cluster to the standard Dirichlet draws does not
aﬀect the analysis and we ignore this for ease of exposition. We will now discuss
our inference algorithm, whose consistency results are presented in Section 3.
Let P = VEVT be the top-K eigendecomposition of P = E[A|Θ]. Typically,
if one assumes B is positive deﬁnite, then ΘB1/2 equals VE1/2 Q for some rank-K
orthogonal matrix Q. Hence, one tries to ﬁrst ﬁnd Q, and then use a post-processing
step to ﬁnd Θ and B. However, in our setting, B is full rank but not necessarily
positive deﬁnite, thereby blocking the avenue. Instead, we proceed from a rather
interesting observation that the population eigenvectors lie on a rotated and scaled
simplex, as shown next.
Lemma 2.3. Let V be the top K eigenvectors of P. Then, under Assumption 2.2,
V = ΘVP , where VP = V(I, :) and I is the indices of rows corresponding to K
pure nodes, one from each community.
Proof. Without loss of generality, reorder the nodes so that Θ(I, :) = I. Then,
VP EVTP = P(I, I) = ρB, so VP ∈ RK×K is full rank. Now, observe that
VP EVT = P(I, :) = ρΘ(I, :)BΘT = ρBΘT . Hence, V = PVE−1 = ρΘBΘT VE−1 =
ΘVP EVT VE−1 = ΘVP .
P

T
Note that Vi := eTi V =
j θij (ej VP ), i.e. the rows Vi of V are convex
combinations of the rows of VP . Thus, the rows Vi lie on a rotated and scaled
simplex, whose corners correspond to the pure nodes. Also, the row norm of Vi is
upper bounded by the maximum norm of the rows in VP , that is, the corners of the
simplex have the highest norm.
This structure immediately suggests the following algorithm. First, ﬁnd the
corners of the rotated and scaled simplex. This can be accomplished using existing
corner-ﬁnding methods such as successive projection algorithm (SPA) [12]. SPA ﬁrst
ﬁnds the node with the maximum norm: i = arg maxj keTj Vk. This node is added
to the set of pure nodes. Then, all remaining rows of V are projected on to the
subspace that is orthogonal to the span of the pure nodes. The process is repeated
for K iterations, and yields a set SP of K pure nodes, one from each community.
Now, we set VP = V(SP , :). Finally, following Lemma 2.3, we can estimate Θ by
T
VV−1
P and B by VP EVP . Since Θ(SP , :) is a permutation matrix, these estimates
are correct up to a permutation.

7

Algorithm 1 Prune
Rn×K ,

Input: Empirical eigenvectors V̂ ∈
an integer r, and two numbers
q, ε ∈ (0, 1).
Output: Set S of nodes to be pruned.
1: for i ∈ n do
2:
vi = keTi V̂k
3: end for
4: S0 = {i : keT
i V̂k ≥ quantile(v, q)}
5: for i ∈ S0 do
6:
di :={Dist. to r nearest neighbors}
P
7:
xi = j dij /r
8: end for
9: S = {i : xi ≥ quantile(x, 1 − ε)}

Algorithm 2 SPACL
Input: Adjacency matrix A, number
of clusters K
Output: Θ̂, B̂, ρ̂.
1: Get the top-K eigen-decomposition
of A as V̂ÊV̂T .
2: S = Prune(V̂, 10, .75, .95)
3: X = V̂([n] \ S, :)
4: Sp = SPA(X T )
5: Xp = X(Sp , :)
6: Θ̂ = V̂Xp−1 .
7: Θ̂ = diag(Θ̂+ 1K )−1 Θ̂+
8: B̂ = Xp ÊXpT
9: ρ̂ = maxi,j B̂ij . B̂ = B̂/ρ̂

In practice, we only observe the empirical eigenvectors, which are noisy versions
of the population eigenvectors. The main idea of pruning is to identify and remove
the nodes which are far away from the population simplex. Note that as the graph
gets sparser, there would be more variance, leading to points having large deviations
from the population simplex. The pruning algorithm (Algorithm 1) essentially ﬁnds
these by ﬁrst ﬁnding contenders of pure nodes, i.e. nodes i whose eigenvector rows
V̂i := eTi V̂ have high norm. Among these, it identiﬁes points which do not have
too many nearest neighbors, or in other words, have larger average distance to their
nearest neighbors in comparison to others.
Our main algorithm, called “Sequential Projection After CLeaning” , i.e. SPACL
(Algorithm 2), ﬁrst prunes away noisy high-norm points, then applies SPA to ﬁnd
corners, and ﬁnally estimate Θ and B.
Population simplex
Empirical simplex
Nodes pruned

0.05

0.75

Without pruning
With pruning

0.7

0.65

0

0.6

0.55

0.5
-0.05
0.45

0.4
0

0.06

0.04

0.02

0

-0.02

-0.04

-0.06

8

(A)

8.5

9

9.5

10

10.5

11

11.5

12

(B)

Figure 1: (A) Nodes picked out by Pruning. (B) Eﬀect of pruning on estimating Θ̂.
8

Here we use a simulated network to show the beneﬁts of pruning. In Figure 1 (A)
we plot the population simplex (red circles), the empirical simplex (green dots) (upto a rotation), and the nodes pruned (blue stars) for a MMSB model with n = 5000,
α = (0.4, 0.4, 0.4), B = (1 − q)I3 + q13 1T3 with q = 0.001 and ρ = 0.007. We can
see that after pruning the nodes remaining are closer to the population simplex and
should lead to more stable estimation. Figure 1 (B) varies ρ from 0.005 to 0.007
leading to average degrees increasing from 8 to 12, and shows the eﬀect of pruning
(blue dashes) over not pruning (red solid line) on the relative estimation error of Θ.

3

Main results

We want to prove that the sample-based estimates Θ̂, B̂ and ρ̂ concentrate around
the corresponding population parameters Θ, B, and ρ. This can be easily done by
applying the theoretical result from [12] which essentially guarantees that the pure
nodes picked using SPA are close to the true pure nodes as long as each of the input
data-points are close to their population counterpart, which lies on a simplex. In our
case, the input data points are the rows V̂i of the matrix V̂ ∈ Rn×K , comprising the
top-K eigenvectors of A. Existing techniques like the Davis-Kahan Theorem [32]
show the convergence of empirical eigenvectors to their (rotated) population counterpart, but they only provide convergence in the Frobenius norm kV − V̂OkF or
the operator norm kVVT − V̂V̂T k. These translate into loose bounds on the rows
of V̂. Otherpexisting techniques [6, 5, 20] can be applied to show that rows of V̂
have ÕP (1/ nρ2 ) relative error, but this is only meaningful when the degree grows
faster than square root of n, i.e. the dense degree regime.
√
We are able to show that keTi (V̂V̂T −VVT )k/keTi (VVT )k = ÕP (1/ nρ) (Corollary 3.3). Thus, the row-wise relative eigenvector deviation converges to zero even
when average degree only grows faster than the logarithm of n. This requires a subtle argument which uses Kato’s contour integration technique [22]), combined with
a derivation of delocalization properties of the population eigenvectors (i.e., each
√
element of the matrix V has order 1/ n), and path counting arguments. While we
prove row-wise bounds on the n × n projection matrix V̂V̂T corresponding to the
empirical eigenvectors, we also show that both Algorithms 1 and 2 return the same
set of pure nodes whether we invoke them with V̂V̂T or V̂ (Lemma 3.4).
It may seem that the rest of the analysis is straightforward via a simple use of
existing results from [12]. These results guarantee that when the input data-points
are close to their population counterparts lying on a possibly rotated and scaled
simplex, the SPA algorithm returns nearly pure nodes or corners. This immediately
also provides convergence guarantees on B̂ and Θ̂.
√
The row-wise eigenvector bounds also enables us to provide Õ(1/ nρ) rates of
convergence on θ̂i . This is analogous to strong consistency in clustering, and to our
knowledge the ﬁrst such result for detecting mixed memberships.
However, a traditional analysis of empirical eigenvectors would yield a subopti9

p

mal ÕP (1/ nρ2 ) rate on these parameter estimates, thereby proving consistency
only in the regime where average degree grows faster than square root of n, not logarithm of n [23]. Furthermore, in a traditional analysis the error bound is inversely
proportional to λ∗ (B)2 , whereas our bounds are only inversely proportional to λ∗ (B).
Recall that, in the simple setting of ρB = (pn −qn )IK +qn 1K 1TK , λ∗ (B) = (pn −qn )/ρ,
i.e. it captures the separation between the connection probabilities within and across
clusters.
We achieve this by a new construction in which we consider groups of population eigenvalues lying within specially constructed intervals, such that the ratio of
the largest and smallest eigenvalues within any interval is controlled. In contrast,
concurrent work [1, 10, 28] either has a relatively worse dependence on λ∗ (B) or
implicitly assumes that the population eigenvalues are of the same order. Note that,
if the population eigenvalues are of the same order, average expected degree in [1]
can be a constant times log n, whereas we require it to grow faster than log2 n.
Theorem 3.1 (Row-wise convergence). Let ν = α0 / mini αi . If Assumption 2.1 is
satisfied, then for a fixed i ∈ [n], with probability at least 1 − O(Kn−3 ),
eTi (V̂V̂T

T

− VV ) = Õ

K 2.5 ν 1.5 (1 + α0 )1.5
√
n ρλ∗ (B)

!

.

Remark 3.1 (Connection to row-wise error). Note that the above row-wise error
immediately gives us an error bound on rows of V̂.
keTi (V̂ − V(VT V̂))k = keTi (V̂V̂T − VVT )V̂k ≤ keTi (V̂V̂T − VVT )k
The K × K matrix VT V̂ takes out the projection of V on V̂ from V̂.



Remark 3.2 (Connection to Stochastic Blockmodels). While we prove this for the
Mixed Membership Stochastic Blockmodel, the same result extends to Stochastic
Blockmodels as well, by taking α0 → 0. We believe it should also be applicable to
other random graph models when the average degree grows faster than the logarithm
of n, as long as the population eigenvectors are delocalized, and the rank of the
expectation matrix grows slowly compared to the size of the matrix.
Typically, in the Spectral Clustering literature, one proves weak consistency [27,
19], or exact recovery via a two step procedure [21]. Our analysis of eigenspaces
immediately shows a path to establish strong consistency without a second step.
We show this in the Stochastic Blockmodel with simple equal size cluster setting
with ρB = (pn − qn )I2 + qn 12 1T2 . Assume without loss of generality that the nodes
are arranged such that the ﬁrst half of the nodes are from the ﬁrst cluster. Here the
√
√
ﬁrst eigenvector is 1n / n and the second eigenvector is [1Tn/2 |− 1Tn/2 ]T / n. We also
have α0 → 0, K = ν = 2, and λ∗ (B) = (pn −qn )/ρ. Our results imply exact recovery
√
√
with probability greater than 1 − O(K/n2 ), as long as (pn − qn )/ ρ = Ω(1/ n) up
to logarithmic factors.

10

One of the elements of the proof of Theorem 3.1 is a bound on eTi (A − P)Vej .
A naive
√ application of Azuma’s concentration inequality shows that this quantity is
OP ( log n). This is because the above is a sum of n independent centered random
variables, whose squared norms sum to at most one, since Vej has unit norm. This
unfortunately is not dependent on the sparsity parameter ρ (same scale as the standard deviation of elements of A). If the j th population eigenvector has a one in one
location and zeros elsewhere, this quantity is in fact O(1). However, one can bring
in the sparsity parameter ρ, if Vej can be shown to be delocalized, i.e. each entry
√
is of the order 1/ n.
Lemma 3.2 (Delocalization of population eigenvectors). We have that, maxi keTi Vk2 ≤
2
T
1/λK (ΘT Θ) and mini keTi VP k2 ≥ 1/λ1 (ΘT Θ). Furthermore, max

 i kei Vk ≤
2ν(1+α0 )
n
T Vk2 ≥ 2 with probability larger than 1−K exp −
and
min
ke
i
i
n
3n
36ν 2 (1+α0 )2
respectively.
The above result immediately gives a high probability upper bound on kVej k∞ .
Furthermore, since keTi VVT k = keTi Vk, this and Theorem 3.1 show that the relative
√
error of row-wise deviation of eigenspaces decays at the ÕP (1/ nρ).
Corollary 3.3 (Row-wise relative convergence). Let ν = α0 / mini αi . If Assumption 2.1 is satisfied, then for a fixed i ∈ [n],
eTi (V̂V̂T − VVT )
keTi VVT k

= Õ

K 2.5 ν 1.5 (1 + α0 )1.5
√
nρλ∗ (B)

!

with probability larger than 1 − O(Kn−3 ).
Lemma 3.2, in conjunction with a Bernstein type argument, shows that eTi (A −
√
ρ log n (Lemma 5.4). Note that this is only one element
P)Vej is in fact OP
of the full analysis, which also involves bounding projections of diﬀerent powers of
A − P on the population eigenspace. We defer it to Section 5.
So far we have talked about row-wise bounds on empirical eigenspaces. But it
seems cumbersome to apply our algorithms on the n×n V̂V̂T matrix. The following
simple result shows that our algorithms return the same set of pure nodes using V̂
and V̂V̂T . Thus, for the algorithm we simply use V̂.
Lemma 3.4. The pruning algorithm (Algorithm 1) and the SPA algorithm will
return the same node indices on both V̂ (V̂T for SPA) and V̂V̂T .
Concentration of (Θ̂, B̂, ρ̂). Once the pure nodes are recovered, SPACL estimates Θ̂ and B̂ using V̂. Now we will show that these estimates are also consistent.
The proof technique here is relatively straightforward. An application of Theorem
3 in [12] shows that the nodes found by Algorithm 2 are almost pure. Applying
this result requires some additional conditions on the model parameters, which are
detailed below.
11

Assumption 3.1. The model parameters ν, α0 , ρ, λ∗ (B), K and the number of
1.5
3 2.5
0)
.
nodes should satisfy λ∗ (B) = Ω̃ K ν √(1+α
nρ
With these added assumptions, analyzing the relative error of Θ̂ and B̂ still
requires careful analysis. This is because since a naive application of existing results
from random matrix theory will show that the error bounds depend on the inverse
of the square of λ∗ (B). Here again the carefully constructed intervals of population
eigenvalues help to make this dependence proportional to just the inverse of λ∗ (B).
Theorem 3.5. Let Θ̂ be obtained from Step 6 of Algorithm 2. If Assumptions 2.1, 2.2,
and 3.1 are satisfied, there exists a permutation matrix Π, such that with probability
at least 1 − O(K/n2 ), ∀i ∈ [n],


eTi Θ̂ − ΘΠ
eTi Θ



= Õ

kρ̂B̂ − ρΠT BΠkF
= Õ
ρkBkF

K 3.5 ν 3 (1 + α0 )2
√
ρnλ∗ (B)
K 3 ν 3 (1 + α0 )2
√
ρnλ∗ (B)

!

!

,

(1)
(2)

Remark 3.3. The above theorem essentially shows that when degree grows faster
than logarithm of n, SPACL consistently estimates the soft community membership
vector for each node with high probability, which is akin to strong consistency in
clustering. In contrast, existing work shows that suitably normalized variants of
kΘ̂ − ΘkF or kΘ̂(:, i) − Θ(:, i)k1 disappears for i ∈ [K], which is analogous to weak
consistency in clustering. We explain the implication of the above results in the
context of a simple model where α = α0 1K /K and ρB = (pn − qn )IK + qn 1K 1TK .
√
If K is Θ(1), the smallest singular value of B can shrink as λ∗ (B) = Ω̃(1/ nρ). In
the traditional form of lower bound on separation probabilities, if we have (pn −
√
√

qn )/ pn = Ω̃(1/ n), the above condition on λ∗ (B) will also be satisﬁed.

4

Experimental results

We present both simulation results and real data experiments to compare SPACL
with existing algorithms for overlapping network models. We compare with the
Stochastic Variational Inference algorithm (SVI) [13], a geometric algorithm for nonnegative matrix factorization for MMSB models with equal Dirichlet parameters
(GeoNMF) [20], Bayesian SNMF (BSNMF) [25], the OCCAM algorithm [33] for
recovering mixed memberships, and the SAAC algorithm [18].1
1

We were not able to run Anandkumar et al.’s GPU implementation [3] because a required library
CULA is no longer open source, and we could not get good results with the CPU implementation
with default settings.

12

(A)

(B)

(C)

(D)

Figure 2: (A) Error against ǫB : β = (0.5 − ǫB , 0.5, 0.5 + ǫB ). We use diag(B) =
β/ max(β) and Bij = 0.05 for i 6= j. (B) Error against ǫB : similar to (A) with
βi = 0.5 + (i − 4)ǫB , i ∈ [K], K = 7, Bij = 0.2 for i 6= j (C) Error against ǫ: we use
Bii = 1, Bij = ǫ for i 6= j and K = 7. (D) Error against λK (B).
For real data experiments we collect two large datasets (with up to 100,000
nodes) from the DBLP corpus and compare SPACL with GeoNMF, SVI, OCCAM,
BSNMF and SAAC. One of these is assortative (B has positive eigenvalues) and one
which is dissortative (B has negative eigenvalues). We show that for the dissortative
setting, SPACL vastly outperforms other methods. For datasets simulated from the
MMSB model, we only compare with SVI and GeoNMF, since the others make
diﬀerent model assumptions.

4.1

Simulations

On simulated data, we compare the relative error (kΘ̂ − ΘkF /kΘkF ) of diﬀerent
methods, averaged over 10 random runs in a range of parameter settings. In particular, we investigate their sensitivity to the cluster-connection matrix B, and ﬁnally
to negative eigenvalues of B. The largest row-wise relative error has similar trends.
By default, we set n = 5000. We include more simulation results involving the
Dirichlet parameter α, the number of communities K in the Appendix (Section E).
Changing B: First, we set K = 3, ρ = 0.2, α = (1/3, 1/3, 1/3), β = (0.5 −
ǫB , 0.5, 0.5 + ǫB ). We use diag(B) = β/ max(β) and Bij = 0.05 for i 6= j. The result
is shown in Figure 2 (A). We see that SPACL and GeoNMF are comparable and
signiﬁcantly better than SVI. In order to make the problem harder, we change K = 7,
ρ = 0.15, αi = 0.1, βi = 0.5 + (i − 4)ǫB , i ∈ [K]. We use diag(B) = β/ max(β) and
Bij = 0.2 for i 6= j. As shown in Figure 2 (B), in this harder regime, SVI does not
converge, and SPACL outperforms GeoNMF. Also GeoNMF has larger variability
than SPACL.
In Figure 2 (C), we plot the relative error in estimating Θ against increasing oﬀ
diagonal noise ǫ of B. We take K = 7, ρ = 0.15, αi = 1/3, Bii = 1, i ∈ [K]. We
see that SVI performs poorly, and while GeoNMF and SPACL have a similar trend,
13

Table 2: Network statistics. Mono corresponds to author-author graphs, whereas
Bi corresponds to bipartite paper-author graphs.
Dataset
# nodes n
# communities K
Average Degree
Overlap %

DBLP1
Mono
Bi
30,566 103,660
6
12
8.9
3.4
18.2
6.3

DBLP2
Mono
Bi
16,817 50,699
3
6
7.6
3.4
14.9
5.6

DBLP3
Mono
Bi
13,315 42,288
3
6
8.5
3.6
21.1
5.7

DBLP4
Mono
Bi
25,481 53,369
3
6
5.2
2.6
14.4
6.9

DBLP5
Mono
Bi
42,351 81,245
4
8
6.8
3.0
18.5
9.7

SPACL outperforms GeoNMF for large ǫ.
Changing λK (B): We conclude the simulations with experiments on B with negative eigenvalues. We generate
B so that the smallest
eigenvalue λK (B) of B is


1
0.2
0.1


0.5
0.075 · i and vary i ∈ [15]. As i grows,
negative. We set B = 0.2
0.1 0.075 · i
0
λK (B) becomes more negative. We set K = 3, ρ = 0.15, α = (1/3, 1/3, 1/3). We
plot the relative error against λK (B) in Figure 2 (D). We see that both GeoNMF
and SVI perform poorly, whereas SPACL performs better by a large margin.

4.2

Real Data

We use the two types of DBLP networks obtained from the DBLP dataset2 , where
each ground truth community is a group of conferences on one topic. The authorauthor networks were constructed in [20], whereas we construct the bipartite networks by using both papers and authors as nodes. Each community is split into two,
the paper community and the author community. The papers are pure nodes since
they belong to one conference and hence one community, whereas the authors may
belong to more than one community, since they often publish in many conferences.
The details of the subﬁelds can be found in [20]. We have two simple preprocessing
steps for the adjacency matrix: 1) delete nodes that do not belong to any community; 2) delete nodes with zero degree. The statistics of the network are in Table
2, which show that despite being sparse, the networks have large overlaps between
communities. The amount of overlap is measured by the number of overlapping
nodes divided by n.
Implementation details: For real world networks, specially the bipartite networks, when average degree of graphs with 100,000 nodes is smaller than four, some
nodes may have extremely small values of Θ̂ and the corresponding rows may in fact
become zero after thresholding. For those we essentially cannot make any prediction.
This is why for Step 7 of Algorithm 2, we threshold all values smaller than 10−12
to zero and we do not normalize rows which are all zeros. This does not make any
diﬀerence for simulations, but for the real world networks, this stabilizes the results.
2

http://dblp.uni-trier.de/xml/

14

SPACL
GeoNMF

SVI
BSNMF

OCCAM
SAAC

0.30

RCavg

0.25
0.20
0.15
0.10
0.05
0.00

DBLP1

DBLP2

DBLP3

DBLP4

10 3
10 2
10 1
10 0

DBLP1

DBLP2

DBLP3

DBLP4

DBLP5

(B) Running time (log scale) on DBLP.

OCCAM
SAAC

0.4

SPACL
GeoNMF

10 5

SVI
BSNMF

OCCAM
SAAC

10 4

0.3

RCavg

10 3

0.2

10 2

0.1
0.0

OCCAM
SAAC

Running time /s

SVI
BSNMF

SVI
BSNMF

10 4

10 -1

DBLP5

(A) Rank correlation on DBLP.
SPACL
GeoNMF

SPACL
GeoNMF

10 5

Running time /s

0.35

10 1

DBLP1

DBLP2

DBLP3

DBLP4

DBLP5

DBLP1

(C) Rank correlation on DBLP
bipartite graphs.

DBLP2

DBLP3

DBLP4

DBLP5

(D) Running time (log scale) on DBLP
bipartite graphs.

Figure 3: Results on DBLP networks. For bipartite graphs, BSNMF was out of
memory for DBLP1 and DBLP5.
Evaluation Metric: For author nodes, we construct the corresponding row of
Θ by normalizing the number of papers an author has in diﬀerent ground truth
communities. We present the averaged Spearman rank correlation coeﬃcients (RC)
between Θ(:, a), a ∈ [K] and Θ̂(:, σ(a)), where σ is a permutation of [K]. The
formal deﬁnition is:
RCavg (Θ̂, Θ) =

K
X
1
max
RC(Θ̂(:, i), Θ(:, σ(i))).
K σ i=1

Note that RCavg (Θ̂, Θ) ∈ [−1, 1], and higher is better. Since SAAC returns
binary assignment, we compute its RCavg against the binary ground truth.
Performance: We report the RCavg score in Figure 3. We see that SPACL has
better prediction accuracy on the bipartite versions of the datasets. This is because
the bipartite author-paper network retains information that is lost when the authorauthor networks are constructed. Also SPACL outperforms all other networks on
bipartite networks, since these are dissortative and the corresponding B will have
negative eigenvalues. On co-authorship graphs, SPACL performs comparably to
15

GeoNMF, while the other methods are worse. Both SPACL and GeoNMF are much
faster than the competing algorithms.

5

Analysis

We ﬁrst present the proofs of the identiﬁability results, followed by the proof of the
row-wise eigenvector bounds.

5.1

Identifiability

Proof of Theorem 2.1. Without loss of generality, we absorb ρ in B, and reorder
nodes so that the ﬁrst K nodes contain one pure node from each community. Thus,
Θ1:K,: = IK .
Let P = VEVT be the eigen-decomposition of P, with V ∈ Rn×rank(B) . Let
VP = V1:K,:. Lemma 2.3 shows that V = ΘVP . Thus, for any node i, Vi lies in
the convex hull of the K rows of VP , that is, Vi ∈ Conv(VP ). We will slightly
abuse the classical notation to denote by Conv(M) the convex hull of the rows of
matrix M.
Now, suppose P can be generated by another set of parameters (Θ′ , B′ ), where
Θ′ has a diﬀerent set of pure nodes, with indices I =
6 1 : K. By the previous
argument, we must have VI ⊆ Conv (VP ). Since (Θ′ , B′ ) and (Θ, B) have the same
probability matrix P, they have the same eigen-decomposition up to a permutation
of the communities. Thus, swapping the roles of Θ and Θ′ and reapplying the
above argument, we ﬁnd that VP ⊆ Conv (VI ). Then Conv (VP ) ⊆ Conv (VI ) ⊆
Conv (VP ), so we must have Conv (VP ) = Conv (VI ). This means the pure nodes in
Θ and Θ′ are aligned up to a permutation, that is, VI = MVP , where M ∈ RK×K
is a permutation matrix.
Now, V = ΘVP = Θ′ VI = Θ′ MVP , which implies (Θ − Θ′ M)VP = 0. Since
V = ΘVP and rank(Θ) = K, we have rank(VP ) = rank(V) = rank(B). Hence, if
rank(B) = K, VP is full rank, so Θ = Θ′ M. Thus, Θ and Θ′ are identical up to a
permutation. To have the same P, B and B′ must also be identical up to the same
permutation. Hence, the MMSB model is identiﬁable. This proves part (a).
Now, suppose rank(B) = K − ℓ < K. We ﬁrst permute the columns of Θ, and
the rows and columns of B, so that
B=

"

C

CW
WT CW

WT C

#

,

where C ∈ R(K−ℓ)×(K−ℓ) is full rank, and W ∈ R(K−ℓ)×ℓ . We see that
C
WT C

h

h

IK−ℓ W
IK−ℓ W

i
i

= V1:(K−ℓ),: EVTP ,
= V(K−ℓ+1):K,: EVTP .

16

(3)

The ﬁrst equation shows that rank(V1:(K−ℓ),: ) = rank(C) = K − ℓ, so V1:(K−ℓ),: is
full rank. Hence,
T

V(K−ℓ+1):K,: = W V1:(K−ℓ),: ⇒ VP =
⇒V = Θ

"

#

IK−ℓ
WT

"

IK−ℓ
WT

#

V1:(K−ℓ),:

V1:(K−ℓ),: .

(4)

Case 1: rank(B) = K − 1 (so W is a vector) and WT 1K−ℓ 6= 1.
We prove identiﬁability by contradiction. Suppose a diﬀerent MMSB model
(Θ′ , B′ ) yields the same probability matrix P, and hence the same eigenvector matrix
V (up to a column permutation). Then,
Θ

"

IK−ℓ
WT

#

′

V1:(K−ℓ),: = V = Θ

"

IK−ℓ
WT

#

V1:(K−ℓ),:

⇒Θ′ = Θ.

(5)

The above equation is derived using Θ1K = Θ′ 1K = 1n , and WT 1K−ℓ 6= 1.
Clearly B′ = B as well, so the MMSB model is identiﬁable. From Equation (3),
we have B(K−ℓ+1):K,: = WT B1:(K−ℓ),: , so WT 1K−ℓ 6= 1 iﬀ the last row of B is not
a aﬃne combination of the remaining rows. It is easy to see that the same holds for
any row of B. This proves part (b).
Case 2: rank(B) = K − 1 and WT 1K−ℓ = 1, or rank(B) < K − 1.
We will construct a Θ′ 6= Θ that yields the same probability matrix P. Let the
mixed node be m, so θmj > 0 for all communities j. We use
θj′

=

(

θj
6 m
h
i if j =
,
T
T
θm + ǫβ −W | Il if j = m,

′
∈ (0, 1) for all communities j, and β ∈ Rℓ 6= 0 is
where ǫ is small enough that θmj

h

i

such that β T −WT 1K−ℓ + 1ℓ = 0. Note that such a β always exists when ℓ > 1

and can be arbitrary vector when WT 1K−ℓ = 1ℓ . Hence, each row of Θ′ sums to 1,
and Θ′ is a valid community-membership matrix. Additionally, Θ′ VP = ΘVP .
Finally, we will show that (Θ′ , B) and (Θ, B) generate the same probability
matrix. Note that B = P1:K,1:K = VP EVTP . Hence,
ΘBΘT = ΘVP EVTP ΘT = VEVT = P = Θ′ VP EVTP Θ′T = Θ′ BΘ′T .

This proves part (c).
Proof of Theorem 2.2. Consider an MMSB model parametrized by (Θ(1) , B(1) ), with
T
P = Θ(1) B(1) Θ(1) (we absorb ρ in B without loss of generality). We want to
17

construct a (Θ(2) , B(2) ) that gives the same probability matrix P. The idea is to
construct a matrix M such that Θ(2) = Θ(1) M and B(2) = M−1 B(1) (MT )−1 . The
diﬃculty is in ensuring that all constraints are satisﬁed: Θ(2) 1K = 1n , Θ(2) ≥ 0,
(2)
and 0 ≤ Bij ≤ 1 for all i, j.
Without loss of generality, suppose that the ﬁrst community does not have any
(1)
pure nodes. In other words, for all nodes i ∈ [n], θi1 ≤ 1 − δ for some δ > 0.
Consider the following M:
M=

"

1 + (K − 1)ǫ2
0

−ǫ2 1TK−1
ǫ1K−1 1TK−1 + (1 − (K − 1)ǫ)IK−1

#

,

where ǫ is a small positive number (0 < ǫ < δ). It is easy to check that M is full
rank (for small enough ǫ) and M · 1K = 1K . Hence, Θ(2) 1K = Θ(1) M1K = 1n
For any node i and for j > 1,
(2)

(1)

θi1 = θi1 (1 + (K − 1)ǫ2 ) ≥ 0,
(2)

(1)

θij = −θi1 ǫ2 +

K
X
ℓ=2

(1)

(1)

θiℓ Mℓj ≥ −θi1 ǫ2 + ǫ

K
X
ℓ=2

(1)

θiℓ ≥ ǫδ2 > 0,

(1)

where we used ǫ < δ and θi 1K = 1. Hence, Θ(2) ≥ 0.
Finally, we must show that B(2) = M−1 B(1) (MT )−1 has all elements between 0
and 1. Note that
M − IK =

"

(K − 1)ǫ2
0

−ǫ2 1TK−1
ǫ1K−1 1TK−1 − (K − 1)ǫIK−1

#

,

so kM − IK kF → 0 as ǫ → 0. Since M−1 is continuous in M, we have M−1 − IK
0. Thus,
B(2) − B(1)

F

= M−1 B(1) (MT )−1 − B(1)
≤ M−1 − IK

2

F

B(1)

F

F

→

F

+ 2 M−1 − IK

F

B(1)

F

→ 0 as ǫ → 0.
(1)

(2)

Since Bij ∈ (0, 1), we have Bij ∈ (0, 1) for ǫ small enough, completing the proof.

5.2

Eigenspace Row-wise Concentration

Before presenting the analysis of the row-wise error-bounds of empirical eigenvectors,
we present a discretization scheme of the population eigenvalues, which later helps
in getting a better dependence of the overall row-wise error on the smallest singular
value of P, which can also be thought of as the separation between blocks.
18

Definition 5.1 (A discretization of eigenvalues). Lets divide the eigenvalues of P
up into the positive ones (S + ) and negative ones (S − ). We will deﬁne the dividing
process only for the positive eigenvalues. The same argument can be carried out on
the negative eigenvalues by taking their absolute values. We start with the smallest
eigenvalue in S + . Denote this by λ∗+ . We use this as gap g1 and keep moving through
the eigenvalues in S + in increasing order until we ﬁnd two consecutive eigenvalues
which have gap g2 > λ∗+ . Thus every pair of consecutive eigenvalues in the kth
interval is within gk distance. Furthermore, as k grows, the gap increases. We
deﬁne sk and ek as the starting and ending index of eigenvalues of the kth interval.
Thus λ∗ (P) ≤ g1 < g2 < g3 . . .
Sk+ = {λsk , . . . λek ∈ S + : λkj − λkj+1 ≤ gk for sk ≤ kj ≤ ek , λsk −1 − λsk > gk }

Deﬁne by VSk the columns of V corresponding to Sk . Note that nk := |Sk | ≤ K,
k nk = K, and the number of intervals (I) can also be at most K. Let the number
of intervals with positive eigenvalues be I + and the ones with negative eigenvalue
be I − . Figure 4 shows an example for S + when there are 10 positive eigenvalues.

P

λ10 λ9

0
g1

λ8

λ6

g2

S1+
λe1

λ7

λs1

λ5

λ4

λ2 λ1

g3

S2+
λe2

λ3

S3+

λs2

λe3

λs3

Figure 4: An illustration of Deﬁnition 5.1.
The above discretization lets us control the ratio of the largest eigenvalue in each
interval and the gap between an interval and the next as follows.
Lemma 5.1. Consider the intervals defined in Definition 5.1. We have: λsk ≤
i=1 ni gk .

Pk

Proof. We prove this by induction. First, note that the smallest positive eigenvalue
is larger than λ∗ (P) by deﬁnition. For k = 1, λs1 − λe1 ≤ (n1 − 1)λe1 , and hence
P
λs1 ≤ n1 λe1 = n1 g1 . Now assume that λsk ≤ ki=1 ni gk . Hence,
λsk+1 ≤ (nk+1 − 1)gk+1 + (λek+1 − λsk ) + λsk = nk+1 gk+1 +

k
X
i=1

ni gk ≤

k+1
X

ni gk+1

i=1

The last step holds since gk < gk+1 .
We also have the following lower bound on λ∗ (P).
Lemma 5.2. Let ν = α0 / mini αi and let λ∗ (P) denote the K th largest singular
value of P. We have λ∗ (P) ≥ ρλ∗ (B) · λK (ΘT Θ) a.s., and


P λ∗ (P) ≥

nρλ∗ (B)
2ν(1 + α0 )





≥ 1 − K exp −
19



n
.
2
36ν (1 + α0 )2

In order to prove our entrywise concentration bound, we will ﬁrst introduce the
notion of matrix resolvents and useful identities on resolvents.
Definition 5.2. A resolvent of a matrix M ∈ Rn×n is deﬁned as GM (z) = (M −
T
P
i (M)
zI)−1 , where z 6∈ {λi (M)}ni=1 . We can also write the resolvent as ni=1 vi (M)v
λi (M)−z ,
where vi (M) is the ith eigenvector of M.
Let us denote by
Ez = diag



λi
z(λi − z)

K !

Mz = VEz VT

,

(6)

i=1

As we see below this matrix is an integral part of the resolvent of the expectation
matrix P.
GP (z) =

n
X
vi vTi
i=1

λi − z

=

K
X

vi vTi

i=1



1
1
+
λi − z z



−

I
I
= Mz −
z
z

(7)

We will use a standard technique to compute eigenspaces of matrices (also used
in [22] Lemma A.2). Consider an interval (a, b) such that no eigenvalue of a symmetric matrix M equals a √
or b. Now √
consider a √
contour C √
in the complex plane
which passes through a + γ −1, a − γ −1, b − γ −1, b + γ −1 in counter clockwise direction, where γ > 0. From the Cauchy integration formula, we know that
1
√

2π −1

I

C

GM (z)dz = −

X

vi (M)vi (M)T

(8)

i:λi (M)∈(a,b)

We consider a sequence of non-overlapping contours Ck , k ∈ [I] (I ≤ K) created
using ak , bk , γk , where kA − Pk < ak < bk , and none of the eigenvalues of A or P
equal ak , bk for k ∈ [I]. Let Vk denote the n × nk matrix with the eigenvectors of P
corresponding to eigenvalues in (ak , bk ). Similarly deﬁne V̂k denote the eigenvectors
of A corresponding to eigenvalues in (ak , bk ).
Hence, using the Cauchy integration formula (8), we have:
Vk VTk − V̂k V̂Tk =

1
√

2π −1

I

Ck

(GA (z) − GP (z)) dz

(9)

Furthermore, it is not hard to check that, ∀x ∈ [n],
eTx



Vk VTk

−

V̂k V̂Tk



=

1
√

2π −1

I

Ck

eTx (GA (z) − GP (z)) dz

We bound the Frobenius norm of the above quantities using Lemma 5.3.

20

(10)

Lemma 5.3. We have:
eTx

I
X

(Vk VTk − V̂k V̂Tk ) ≤

k=1

where

I
X
bk − ak + 2γk

π

k=1

max(P1 (z) + P2 (z)), (11)
z∈Ck

P1 (z) = |z|kGA (z)kkA − PkkEz kkeTx GA−P (z)Vk,

P2 (z) = keTx GA−P (z)(A − P)VkF kEz k.

Proof. Since GA (z) − GP (z) = GP (z)(P − A)GA (z),




I
GA (z) − GP (z) = Mz −
(P − A)(GA (z) − GP (z)) + GP (z)(P − A)GP (z)
z
Bringing A−P
z (GA (z)− GP (z)) to the LHS, and using the deﬁnition of the resolvent
of A − P we get:







}

GA (z) − GP (z) = zGA−P (z) Mz (A − P)(GA (z) − GP (z)) + GP (z)(A − P)GP (z)
{z

|

}

R0

{z

|

R

(12)

As it turns out, each of the rows of zR either have small Frobenius norm or they
disappear when combined with zGA−P (z) post integration. We will show this step
by step. Note that, using Equation (7), R in the above equation can be decomposed
as:






I
I
(A − P) Mz −
GP (z)(A − P)GP (z) = Mz −
z
z
A−P
A−P
+ Mz (A − P)GP (z) −
Mz
=
2
|
{z
}
z
}
| {z }
| z {z
R1

R2



R3

Next, we show that zGA−P (z)R1 disappears upon integration. But ﬁrst, we
note that since by construction ∀z ∈ Ck , ∀k, |z| ≥ ak > kA − Pk, we have the
following series expansion for GA−P (z),
GA−P (z) = −



1X A−P
z t≥0
z

t

(13)

Since none of the contours contain zero, Equation (13) immediately gives:
1
√

2π −1

I

Ck

zGA−P (z)R1 dz = −

21

XI

t≥1 Ck

1
z



A−P
z

t

dz = 0

(14)

Thus Equations (10), (12) and (14) give us,
1
√

eTx (Vk VTk − V̂k V̂Tk ) = −

I

eT GA−P (z)z(R0 + R2 − R3 )dz
2π −1 Ck x
bk − ak + 2γk
max keTx GA−P (z)z(R0 + R2 − R3 )k
keTx (Vk VTk − V̂k V̂Tk )k ≤
z∈Ck
π
Now we bound each part individually.
keTx GA−P (z)z(R0 + R2 )k = keTx GA−P (z)zMz (A − P)GA (z)k

≤ keTx GA−P (z)zMz kF kA − PkkGA (z)k

(i)

≤ |z|kGA (z)kkA − PkkEz kkeTx GA−P (z)Vk =: P1 (z)

Step (i) uses Equation (6). Finally we also have:
keTx zGA−P (z)R3 k ≤ keTx GA−P (z)(A − P)VkkEz k =: P2 (z)
Thus, the statement of the lemma follows.
Now we need to
1. Deﬁne the contours and events so that the LHS of Equation (11) can cover
the whole eigenspace.
2. Under these events, we would bound P1 (z) and P2 (z) over each contour. For
this we need to bound keTx GA−P (z)vi k, keTx GA−P (z)(A − P)vi k, where vi
denotes the ith column of V. We also need kEz k, kGA (z)k, kGA−P (z)k,
etc. This requires us to bound |eTi Ht vi | for t ≤ log n, where H := (A −
√
P)/ nρ. For t = 1, we prove the following lemma, which uses the fact that
V is delocalized with high probability (see Lemma 3.2).
Lemma 5.4. Let vk denote the kth population eigenvector of P. If Assumption 2.1
is satisfied,




P ∃k ∈ [K], |eTi Hvk | ≥ 4 log nkvk k∞ = O(K/n3 )
Proof. For ease of notation we will ﬁrst prove this for one population eigenvector
√
v. Recall that H := (A − P)/ nρ. Let Xj = (Aij − Pij )vj , where vj is the j th
component of v. We have |Xj | ≤ kvk∞ =: M . Conditioning on P, i.e. conditioning
on Θ (since B is a ﬁxed parameter) Xj are mean zero independent random variables.
Also note that, since kvk = 1 and Pij ≤ ρ,
X
j

var(Xj |Θ) =

X
j

h

i

E (Aij − Pij )2 vj2 =
22

X
j

Pij (1 − Pij )vj2 ≤ ρ

An application of Bernstein’s inequality gives us:


P |

X
j



t2
Xj | ≥ t Θ ≤ 2 exp − P
2( j var(Xj |Θ) + tM/3)

!

=: 2 exp(−A)

First note that the RHS of the above equation is a decreasing function of t. We set
√
t = 4 max(M, ρ) log n. Consider the following two cases:
√
Case 1: M > ρ: We have t = 4M log n. Hence,
16M 2 log2 n
exp(−A) ≤ exp −
2ρ + 8/3M 2 log n
Case 2: M ≤

!

16M 2 log2 n
≤ exp −
2M 2 + 8/3M 2 log n

!

≤

1
n3

√
√
ρ: We have t = 4 ρ log n. Hence,

16ρ log2 n
exp(−A) ≤ exp −
√
2ρ + 8/3M ρ log n

!

16ρ log2 n
≤ exp −
2ρ + 8/3ρ log n

!

≤

1
n3

Applying this to all K population eigenvectors we have:


2K
√
P ∃k ∈ [K], |eTi (A − P)vk | ≥ 4 max(kvk k∞ , ρ) log n Θ ≤ 3
n

Recall from Lemma 3.2 that, ∀k ∈ [K], kvk k∞ ≤
deﬁned as Ω :=



Θ:

q

1
λK (ΘT Θ)

≤

q

2ν(1+α0 )
n

see that, conditioned on Ω, kvk k∞ ≤
q

√



q

1
.
λK (ΘT Θ)

(15)

Consider the event Ω

. If Assumption 2.1 is satisﬁed, we

ρ, ∀k ∈ [K]. Now we consider the event

A := {∃k ∈ [K], |eTi (A − P)vk | ≥ 4 ρ log2 n}.
P(A) =

Z

Θ∈Ω

P (A|Θ) P(Θ)dΘ +

Z

P (A|Θ) P(Θ)dΘ

Θ6∈Ω



(i) 2K
n
2K
≤ 3 + 1 − P(Θ ∈ Ω) ≤ 3 + K exp −
2
n
n
36ν (1 + α0 )2

 (ii)

≤

3K
n3

Steps (i) and (ii) follow from Lemma 3.2 and the condition on ν from Assumption 2.1
√
respectively. Since ∀k ∈ [K], kvk k∞ ≥ 1/ n, we have:






s

P ∃k ∈ [K], |eTi Hvk | ≥ 4 log nkvk k∞ ≤ P ∃k ∈ [K], |eTi Hvk | ≥ 4


K
=O
n3



2



log n 
n

For 1 < t ≤ log n, we adapt a crucial result from [11], details of which are
deferred to the Appendix (Section C).
23

√
Lemma 5.5. Let H := (A − P)/ nρ. As long as Assumption 2.1 is satisfied for
some constant ξ, for any fixed vector v, for a fixed i and for 1 < t ≤ log n,




P |eTi Ht v| ≤ (log n)tξ kvk∞ |Θ ≥ 1 − exp(−(log n)ξ /3)
We will now deﬁne some events, which will be used extensively to show that
the contours cover all population and empirical eigenvalues, and to bound P1 (z)
and P2 (z) in Equation (11). Let vk be the kth population eigenvector. Under
Assumption 2.1, for t ≤ log n,




n

′
T


E := λK (Θ Θ) ≥


2ν(1 + α0 )





√

 ′′



(i)

P(E¯′ ) ≤ K exp −
(iii)

E := {kA − Pk ≤ C nρ}

n
o


T


e
Hv
E
:=
≤
4
log
nkv
k
,
∀k
∈
[K]

1
∞
k
k
i




n
o


 E := eT Ht v ≤ (log n)tξ kv k , ∀k ∈ [K]
t
∞
k
k
i

P(E¯′′ ) ≤ n−3
(iv)



P(E¯1 ) ≤ O K/n3
(v)

n
2
36ν (1 + α0 )2

 (ii)

≤ K/n3



P(E¯t ) ≤ K exp(−(log n)ξ /3), 1 < t ≤ log n
(16)
In Equation (16) step (i) is proven in the Appendix and step (ii) follows from
Assumption 2.1. For any community membership matrix Θ, P (E¯′′ |Θ) can be
bounded directly using Theorem 5.2 of [19], since Assumption 2.1 guarantees that
nρ = Ω(log n). Hence step (iii) follows. Since P (E¯1 |Θ) and P (E¯t |Θ), 1 < t ≤ log n
can be bounded for all community membership matrices Θ using Lemma 5.4 and
Lemma 5.5 respectively, steps (iv) and (v) follow. To denote order notation condiE′

tioned on event E ′ , we will use, X = O(.) to denote, P(X = O(.)) = P(E ′ ).
Picking the contours Ck : Consider the discretization in Deﬁnition 5.1. For the
kth interval, use γk = gk /4, ak = max(λek − gk /2, (1 + c)kA − Pk), for some c > 0
and bk = λsk + gk /2. If bk ≤ ak , we ignore the contour. If either ak or bk equal an
eigenvalue of A or P, for any ǫ > 0, they can be perturbed by at most ǫ to guarantee
that they do not coincide with eigenvalues of A or P. This is possible because for
a given n, the set {A | A ∈ {0, 1}n×n } is ﬁnite.
Now we bound the spectral norms of GA (z), GP (z), Ez and GA−P (z). Since
the gap between the smallest eigenvalue (in magnitude) of the kth interval and the
largest eigenvalue in the (k − 1)th interval is gk , and by construction (Deﬁnition 5.1)
λ∗ (P) ≤ g1 < g2 < . . . , and λek ≥ gk , we note that for each contour Ck , |z| can be
upper and lower bounded as follows.
|z| ≤

q

b2k + γk2 ≤ bk + γk = λsk + 3gk /4

|z| ≥ max((1 + c)kA − Pk, |λek − gk /2|) ≥ |λek − gk /2| ≥ gk /2

(i)
√
|z − λ̂i | ≥ gk /2 − O( nρ)
 
1
1
1
kMz k = kEz k ≤ max
=O
+
i
λi − z z
gk

|z − λi | ≥ gk /2,

24

(17)
(18)
(19)
(20)

Conditioned on E ′ ∩ E ′′ , ∀i ∈ [n] and ∀z ∈ Ck , Equation (20) follows from Equa√
E ′′
tions (6), (18) and (19). Step (i) in Equation (19), uses |λ̂i − λi | = O( nρ) via
Weyl’s inequality. Finally using Equations (18), (19) and (20) we also have ∀z ∈ Ck ,
conditioned on E ′ ∩ E ′′ ,


1
kGP (z)k = O
gk



X v̂i v̂T
i

kGA (z)k ≤

i

λ̂i − z

1
√
gk − O( nρ)

=O

Now conditioned on E ′ ∩ E ′′ , Equation (21) gives:
kGA (z) − GP (z)k ≤ kGP (z)kkP − AkkGA (z)k = O

√



nρ
O
gk

!

(21)

!

1
√
gk − O( nρ)
(22)

Now we will bound the RHS of Equation (11) in Lemma 5.3.
Lemma 5.6. Let vi denote the ith column of V. Let Assumption 2.1 be satisfied for
some constants c0 and ξ. Consider the events defined in Equation (16). Conditioned
n
′
′′
on ∩log
t=1 Et ∩ E ∩ E ,
|eTx GA−P (z)vi | =



O kvi k∞ + n−2ξ
√

|eTx GA−P (z)(A − P)vi | = O 

λ∗ (P)




.

nρ (log n)ξ kvi k∞ + n−2ξ
λ∗ (P)



.

Proof. We use Equation (13) and Lemmas 5.4 and 5.5. First note that, for H deﬁned
in Lemma 5.5, for 1 ≤ t ≤ log n, conditioned on Et , t ≥ 1,
eTx (A

−
zt

P)t vi



√
nρ(log n)ξ

√
t

(
nρ)
|z|
T t
= ex H vi
≤ 
t

zt
kA−Pk

|z|

t

kvi k∞

t ≤ log n

(23)

t > log n

It is easy to verify that the above holds for t = 0. Conditioned on E ′ ∩ E ′′ , Equation (19), Lemma 5.2 and Assumption 2.1 give:
E′

√
λ (P) ≥ 4 nρ(log n)ξ
∗

√
nρ(log n)ξ E ′ 1
≤
max
k,z∈Ck
|z|
2

⇒

25

(24)

n
Conditioned on E ′ ∩ E ′′ ∩ (∩log
t=1 Et ), Equations (13) and (23) give:

max |eTx GA−P (z)vi | ≤ max

k,z∈Ck

k,z∈Ck

∞ T
ex (A − P)t
1 X
vi
|z| t=0
zt

log n
1 X eTx (A − P)t vi
1 X eTx (A − P)t vi
+
max
≤ max
k,z∈Ck |z|
k,z∈Ck |z|
zt
zt
t=0
t>log n

Eqs (23) and (24)

Eq (24)

kvi k∞
(kA − Pk/|z|)log n+1
+
max
√
k,z∈Ck |z| − nρ(log n)ξ
k,z∈Ck
|z| − kA − Pk
!
√
(2C nρ/λ∗ (P))log n+1
kvi k∞
=O
+
√
√
λ∗ (P)/2 − nρ(log n)ξ
λ∗ (P)/2 − C nρ

≤ max

=



O kvi k∞ + n−2ξ
λ∗ (P)

We also have, for large enough n,






√ 
2C nρ log n+1
λ∗ (P)

≤



C
2(log n)ξ

log n+1

≤ exp(O(log n)−

ξ(log n + 1) log log n) = O 1/(n2ξ ) . Furthermore, using the same argument as before,
max |eTx GA−P (z)(A − P)vi | = max

k,z∈Ck

k,z∈Ck

√

= O

∞ T
X
ex (A − P)t
t=1



zt

vi

nρ (log n)ξ kvi k∞ + n−2ξ
λ∗ (P)

Now we are ready to ﬁnish the proof of Theorem 3.1.




Proof of Theorem 3.1. Our goal is to bound the row norm of VVT − V̂V̂T using
Lemma 5.3. The ﬁrst step is to show:
E ′ ∩E ′′

keTx (VVT − V̂V̂T )k =

I
X

k=1

eTx (Vk VTk − V̂k V̂Tk ) .

(25)

Recall that ak = max(λek − gk /2, (1 + c)kA − Pk). Conditioned on E ′ ∩ E ′′ , Assumption 2.1, and Lemma 5.2, λek − gk /2 ≥ λ∗ (P)/2 = ω(kA − Pk) and hence
ak = λek − gk /2. Hence the intervals are mutually exclusive and cover all the population eigenvalues, proving Equation (25).
n
By triangle inequality, conditioned on E ′ ∩ E ′′ ∩ (∩log
t=1 Et ), from Lemma 5.3, we

26

have:
keTx (VVT
I
X

T

− V̂V̂ )k ≤



I
X

k=1

λsk − λek + 2gk
=
O
gk
k=1

(i)

(ii)

=

I
X

k=1

keTx (Vk VTk − V̂k V̂Tk )k



√

max O
z∈Ck

!

nρ(bk + γk )
keTx GA−P (z)Vk
gk


+keTx GA−P (z)(A − P)Vk
!

!

k
√ X
nρ( ni + 1) keTx GA−P (z)Vk + keTx GA−P (z)(A − P)Vk
O(nk ) max O
z∈Ck



i=1



√
= O(K) O( nρK) max keTx GA−P (z)Vk + max keTx GA−P (z)(A − P)Vk
(iii)

= O

(iv)

= Õ

√
K 3/2 nρ
λ∗ (P)

k,z∈Ck

!

k,z∈Ck

(K + (log n)ξ ) max kvi k∞ + (K + 1)n−2ξ
i

!
1.5



K 2.5 ν 1.5 (1 + α0 )
√
n ρλ∗ (B)

Step (i) uses Equation (17). Step (ii) holds because from construction of the eigenvalue intervals, we have: λsk − λek + 2gk ≤ (nk + 1)gk . It also uses the fact that
P
λsk ≤ gk ki=1 ni by Lemma5.1. Step (iii) follows from Lemma 5.6. Step (iv) uses
maxi kvi k∞ = O

q

ν(1+α0 )
n

= Ω(n−2ξ ). To bound the failure probability, for some

constant ξ > 1 and large enough n, Equation (16) gives:
n
¯′
¯′′
P (E ′ ∩ E ′′ ∩ (∩log
t=1 Et )) ≥ 1 − P (E ) − P (E ) −

6

log
Xn
t=1

P (E¯t ) ≥ 1 − O(Kn−3 ).

Conclusions

This paper proposed a fast and provably consistent algorithm called SPACL for
inferring community memberships of nodes in a network generated by a Mixed
Membership Stochastic Blockmodel (MMSB). The proof of consistency requires a
diﬀerent approach than those used for Stochastic Blockmodels. Our proof has several
new aspects, including a sharp row-wise eigenvector bound using complex contour
integration, a proof of delocalization of the top eigenvectors of the population probability matrix P, and tight bounds on the product of the centered random matrix
A − P with the top eigenvectors of P. We only assume that each community has
at least one pure node, and the community-connection matrix B is full-rank; this is
weaker than the assumption of a positive-deﬁnite B in some prior work. We are also
the ﬁrst to provide both necessary and suﬃcient conditions for identiﬁability under
MMSB, and convergence of the membership vector of each node to its population
27

counterpart. In contrast existing methods establish suﬃcient conditions only and
provide consistency of community memberships as a whole. Results on simulated
and real-world networks shows that SPACL is more accurate and has lower variance
than several related methods, and is fast and scalable as well.

References
[1] E. Abbe, J. Fan, K. Wang, and Y. Zhong. Entrywise eigenvector analysis of random
matrices with low expected rank. arXiv preprint arXiv:1709.09565, 2017.
[2] E. M. Airoldi, D. M. Blei, S. E. Fienberg, and E. P. Xing. Mixed membership stochastic
blockmodels. JMLR, 9:1981–2014, 2008.
[3] A. Anandkumar, R. Ge, D. Hsu, and S. M. Kakade. A tensor approach to learning
mixed membership community models. JMLR, 15(1):2239–2312, 2014.
[4] S. Arora, R. Ge, R. Kannan, and A. Moitra. Computing a nonnegative matrix
factorization–provably. In STOC, pages 145–162. ACM, 2012.
[5] A. Athreya, C. E. Priebe, M. Tang, V. Lyzinski, D. J. Marchette, and D. L. Sussman.
A limit theorem for scaled eigenvectors of random dot product graphs. Sankhya A, 78
(1):1–18, 2016.
[6] S. Balakrishnan, M. Xu, A. Krishnamurthy, and A. Singh. Noise thresholds for spectral
clustering. In NIPS, pages 954–962. 2011.
[7] B. Ball, B. Karrer, and M. E. Newman. Eﬃcient and principled method for detecting
communities in networks. Physical Review E, 84(3):036103, 2011.
[8] J. L. Bentley. Multidimensional binary search trees used for associative searching.
Communications of the ACM, 18(9):509–517, 1975.
[9] A. Beygelzimer, S. Kakade, and J. Langford. Cover trees for nearest neighbor. In
International Conference on Machine Learning, pages 97–104, 2006.
[10] J. Eldridge, M. Belkin, and Y. Wang. Unperturbed: spectral analysis beyond daviskahan. arXiv preprint arXiv:1706.06516, 2017.
[11] L. Erdős, A. Knowles, H.-T. Yau, J. Yin, et al. Spectral statistics of erdős–rényi graphs
i: local semicircle law. The Annals of Probability, 41(3B):2279–2375, 2013.
[12] N. Gillis and S. A. Vavasis. Fast and robust recursive algorithmsfor separable nonnegative matrix factorization. PAMI, 36(4):698–714, 2014.
[13] P. K. Gopalan and D. M. Blei. Eﬃcient discovery of overlapping communities in massive
networks. PNAS, 110(36):14534–14539, 2013.
[14] P. W. Holland, K. B. Laskey, and S. Leinhardt. Stochastic blockmodels: First steps.
Social networks, 5(2):109–137, June 1983. ISSN 0378-8733.
[15] S. B. Hopkins and D. Steurer. Bayesian estimation from few samples: community
detection and related problems. In FOCS, pages 379–390. IEEE, 2017.
[16] J. Jin, Z. T. Ke, and S. Luo. Estimating network memberships by simplex vertex
hunting. arXiv preprint arXiv:1708.07852, 2017.
[17] T. Kato. Perturbation theory for linear operators, volume 132. Springer Science &
Business Media, 2013.
[18] E. Kaufmann, T. Bonald, and M. Lelarge. A spectral algorithm with additive clustering
for the recovery of overlapping communities in networks. In International Conference
on Algorithmic Learning Theory, pages 355–370. Springer, 2016.
[19] J. Lei, A. Rinaldo, et al. Consistency of spectral clustering in stochastic block models.
The Annals of Statistics, 43(1):215–237, 2015.

28

[20] X. Mao, P. Sarkar, and D. Chakrabarti. On mixed memberships and symmetric nonnegative matrix factorizations. In ICML, pages 2324–2333, 2017.
[21] F. McSherry. Spectral partitioning of random graphs. In FOCS, pages 529–537, 2001.
[22] R. I. Oliveira. Concentration of the adjacency matrix and of the laplacian in random
graphs with independent edges. arXiv preprint arXiv:0911.0600, 2009.
[23] M. Panov, K. Slavnov, and R. Ushakov. Consistent parameter estimation in general
stochastic block models with overlaps. arXiv preprint arXiv:1707.01350, 2017.
[24] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes
in C. CUP, CUPA, 2nd edition, 1992.
[25] I. Psorakis, S. Roberts, M. Ebden, and B. Sheldon. Overlapping community detection
using bayesian non-negative matrix factorization. Physical Review E, 83(6):066114,
2011.
[26] A. Ray, J. Ghaderi, S. Sanghavi, and S. Shakkottai. Overlap graph clustering via
successive removal. In 52nd Annual Allerton Conference, pages 278–285. IEEE, 2014.
[27] K. Rohe, S. Chatterjee, and B. Yu. Spectral clustering and the high-dimensional
stochastic blockmodel. The Annals of Statistics, pages 1878–1915, 2011.
[28] M. Tang, J. Cape, and C. E. Priebe. Asymptotically eﬃcient estimators for stochastic
blockmodels: the naive mle, the rank-constrained mle, and the spectral. arXiv preprint
arXiv:1710.10936, 2017.
[29] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of
computational mathematics, 12(4):389–434, 2012.
[30] F. Wang, T. Li, X. Wang, S. Zhu, and C. Ding. Community discovery using nonnegative
matrix factorization. Data Mining and Knowl. Disc., 22(3):493–521, 2011.
[31] X. Wang, X. Cao, D. Jin, Y. Cao, and D. He. The (un) supervised nmf methods for
discovering overlapping communities as well as hubs and outliers in networks. Physica
A: Statistical Mechanics and its Applications, 446:22–34, 2016.
[32] Y. Yu, T. Wang, and R. J. Samworth. A useful variant of the davis–kahan theorem for
statisticians. Biometrika, 102(2):315–323, 2015.
[33] Y. Zhang, E. Levina, and J. Zhu. Detecting overlapping communities in networks using
spectral methods. arXiv preprint arXiv:1412.3432, 2014.

29

Appendix
A

Some Auxiliary Results

Definition A.1. (A construction of rotation matrix) Consider the discretization deﬁned in
Deﬁnition 5.1. The Davis Kahan theorem states that there exists a rotation matrix Ô such
that kV̂−VÔkF is small. In this deﬁnition we will carefully construct this matrix. Consider
the intervals resulting from the discretization of population eigenvalues in Deﬁnition 5.1.
Now, from Theorem 2 of [32], ∃Ôk such that
√
8nk kA − Pk
kRk kF = kV̂Sk − VSk Ôk kF ≤
(26)
gk
Typically the denominator is fk := min(λsk − λsk −1 , min(λek − λek +1 , λek )). We now construct our Ô by stacking the Ôk matrices on the diagonal of a K × K matrix. This is also
a valid rotation matrix. Now, let Ek by the submatrix of E corresponding to eigenvalues in
Sk . Similarly deﬁne Êk . Furthermore, let R := [R1 | . . . |RI ].
Lemma A.1. If Assumption 2.1 holds, then there exists an orthogonal matrix Ô ∈ RK×K
constructed using Definition A.1 that satisfies
√
8KkA − Pk
(27)
kRkF ≤
λ∗ (P)
√ 
Ê − ÔT EÔ = OP K 2 nρ
(28)
F

with probability larger than 1 − O(Kn−3 ).

Proof. Consider the rotation matrix Ô, the residual matrix R constructed as in Deﬁnition A.1. This gives us the Frobenius norm of R as follows, since by construction gk ≥ λ∗ (P).
√
sX
8KkA − Pk
kRk k2F ≤
kRkF ≤
λ∗ (P)
k

Pk
Finally note that, using Lemma 5.1, since λsk ≤ i=1 nk gk ≤ Kgk ,
√
√ √
8nk kA − Pkλsk
T
= OP (K nk nρ)
kRk Ôk Ek kF ≤ kRk kF kEk k ≤
gk
Now we use these intervals as follows.
Ê − ÔT EÔ

F

= V̂ÊV̂T − V̂ÔT EÔV̂T

F

= kAK − (V + RÔT )E(V + RÔT )T kF

≤ kAK − PkF + 2 kRÔT EVT kF + kRÔT EÔRT kF
{z
} |
{z
}
|
P1

p
= OP ( Knρ) + P1 + P2

30

P2

(29)

√
√
KkAK −Pk ≤ K(kA−AK k+kA−Pk) ≤
F ≤
√ The last step is true because kAK −Pk−r
2q Knρ with probability at least 1 − n . As for P1 , note that: P1 ≤ kRÔT EkF ≤
P
3/2 √
T
2
nρ).
k kRk Ôk Ek kF =: OP (K
As for P2 , we have:
 2 
K nρ
|P2 | ≤ kRÔT EkF kRkF = OP
λ∗ (P)
√ √
√
√
Thus the ﬁnal bound is OP (K 2 nρ(max(1/K 3/2 , 1/ K, nρ/λ∗ (P))) = OP (K 2 nρ) using Assumption 2.1. The error probability comes from
the failure of two events, namely

√
ρnλ∗ (B)
∗
kA − Pk = OP ( nρ), and λ (P) = ΩP
. Taking r = 3 we get the required
ν
bound.
Lemma A.2. Let αmax = maxa αa and αmin = mina αa . If ν := α0 /αmin ,




n
n
P λK (ΘT Θ) ≥
≥ 1 − K exp −
2ν(1 + α0 )
36ν 2 (1 + α0 )2




n
2ν(1 + α0 )
≥ 1 − K exp −
P λ1 (VP VTP ) ≤
n
36ν 2 (1 + α0 )2




2α0 (1 + α0 )
n
P λK (VP VTP ) ≥
≥
1
−
K
exp
−
3n (αmax + kαk2 )
36ν 2 (1 + α0 )2


s


n
αmax + kαk2 

≥ 1 − 2K exp −
P κ(VP ) ≤ 3
αmin
36ν 2 (1 + α0 )2

where κ(.) is the condition number of matrix (.).

Proof. Here we consider θi as a random variable. Denote
M̂ = ΘT Θ =

n
X

θi θiT ,

i=1

P

Note that M̂ − E[M̂] =
random matrices. We have
E[θi θiT ] =

i

Xi where Xi are independent mean zero symmetric K × K

diag(α) + ααT
α0 (1 + α0 )

Cov(θi ) =

α0 diag(α) − ααT
.
α20 (1 + α0 )

P
Furthermore, since kθi k1 = 1, and α0 = i αi , we have kXi k ≤ θiT θi + kE[θi θiT ]k ≤
+kαk2
≤ 2. Finally, since the operator norm is convex, Jensen’s inequality gives:
1 + ααmax
0 (1+α0 )
kE[X2i ]k ≤ E[kX2i k] ≤ E[kXi k2 ] ≤ 4. Using standard Matrix Bernstein type concentration
results (Theorem 1.4 of [29]), for large n we get:


t2
P(kM̂ − E[M̂]k ≥ t) ≤ K exp −
=: δt
8n + 4t/3
Now Weyl’s inequality gives, with probability at least 1 − δt ,
|λ1 (M̂) − λ1 (E[M̂])| ≤ t

|λK (M̂) − λK (E[M̂])| ≤ t

31

For the population quantities,
λ1 (E[M̂]) ≤

αmax + kαk2
n,
α0 (1 + α0 )

λK (E[M̂]) ≥

αmin
n
α0 (1 + α0 )


2
n αmax +kαk
n
n
n
2 α0 (1+α0 ) ∈ [ 2ν(1+α0 ) , 2 ] and hence δt ≤ K exp − 36ν 2 (1+α0 )2
αmin
n
n
n
n
2 α0 (1+α0 ) ∈ [ 2ν(1+α0 ) , 2 ]. Hence δt ≤ K exp(− 36ν 2 (1+α0 )2 ).

For λ̂1 (M̂) we take t =
For λ̂K (M̂), we take t =



.

Hence the condition number of M̂ can also be bounded as:
λ1 (ΘT Θ)
≤
κ(Θ Θ) =
λK (ΘT Θ)
T

3n αmax +kαk
2 α0 (α0 +1)
αmin
n
2 α0 (α0 +1)

2

=3

αmax + kαk2
αmin

From Lemma 2.3, V = ΘVP , so
I = VT V = VTP ΘT ΘVP ,
VP is full rank, so ΘT Θ = VP VTP

−1

λK (VP VTP ) =
Similarly we can get λ1 (VP VTP ) ≤

, which gives

1
2α0 (1 + α0 )
≥
.
λ1 (ΘT Θ)
3n (αmax + kαk2 )

2ν(1+α0 )
.
n

And κ(VP ) =

q
p
+kαk2
.
κ(ΘT Θ) ≤ 3 αmaxαmin

From this proof of Lemma 5.2 follows immediately.
Proof of Lemma 5.2. First note that,
q
λ∗ (P) = λ∗ (ρΘBΘT ) = ρλ∗ (BΘT Θ) = ρ λK (ΘT ΘBT BΘT Θ)
q
≥ ρ · λK (ΘT Θ) · λK (BT B) = ρλ∗ (B) · λK (ΘT Θ).

The inequality holds because
∀ full rank symmetric
positive deﬁnite matrix M1 , M2 ∈


−1
,
and
as
λ1 (M−1
M
λ
RK×K , λ1 (M1 M2 )−1 ≤ λ1 M−1
1
1 ) = 1/λK (M1 ) (same for
2
1
M1 and M1 M2 ), we have:
λK (M1 M2 ) ≥ λK (M1 ) λK (M2 ) .
Now, from Lemma A.2, the result follows.
−1
Proof of Lemma 3.2. Note that ΘTΘ = VP VTP  , thus for pure nodes, with Lemma A.2,
n
with probability at least 1 − K exp − 36ν 2 (1+α
,
2
0)
max eTi VP
i

2

= max eTi VP VTP ei ≤ max xT VP VTP x = λ1 VP VTP
i

kxk=1

2ν(α0 + 1)
1
≤
.
=
T
λK (Θ Θ)
n

32



As for other nodes, their rows are convex combinations of the rows of pure nodes and would
be smaller than or equal to the norm of the pure nodes. Thus the result follows.
Now

2
min eTi VP = min eTi VP VTP ei ≥ min xT VP VTP x = λK VP VTP
i

i

kxk=1

1
2α0 (1 + α0 )
2
≥
≥
.
λ1 (ΘT Θ)
3n (αmax + kαk2 )
3n
P
As for node i ∈ [n] that may not be a pure node, as Vi = j θij Vj (W.L.O.G., we make
the ﬁrst K rows of V pure nodes of diﬀerent K communities here), using Jensen’s inequality
we have
=

k

K
X
j=1

θij Vj k2 ≥

K
X
j=1

θij kVj k2 ≥

K
X

θij min eTl VP

j=1

2

l

= min eTl VP

2

l

.

Thus the result follows.

B

Consistency of inferred parameters

Proof of Lemma 3.4. To see that the pruning algorithm returns identical nodes (up-to ties)
is straightforward. This is because the pruning algorithm proceeds by calculating Euclidean
distances between pairs of nodes for nearest neighbor computation. We have
kV̂V̂T (ei − ej )k2 = (ei − ej )T V̂V̂T (ei − ej ) = kV̂T (ei − ej )k2
Thus the pairwise distances between columns of V̂V̂T are the same as that between columns
of V̂T . As for the SPA algorithm, we prove the claim by induction.
Base case: For step k = 1, as
kV̂V̂T ei k2 = eTi V̂V̂T V̂V̂T ei = eTi V̂V̂T ei = kV̂T ei k2 ,
picking max norm will give the same index, denoted as k1 .
Now for V̂T , the vector whose projection is removed is V̂T ek1 , and the normalized vector
is u = V̂T ek1 /kV̂T ek1 k, then for V̂V̂T , the vector whose projection is removed is V̂V̂T ek1
and its normalized vector is u1 = V̂V̂T ek1 /kV̂V̂T ek1 k = V̂V̂T ek1 /kV̂T ek1 k = V̂u.
Now
k(I − u1 uT1 )V̂V̂T ei k2 = k(I − V̂uuT VT )V̂V̂T ei k2

= kV̂(I − uuT )V̂T ei k2 = k(I − uuT )V̂T ei k2 ,

then for step k = 2, picking max norm will also give the same index.
Induction: Suppose for ﬁrst k − 1 ∈ [K − 1] steps SPA on V̂T and on V̂V̂T will
give the same indices as Sk−1 , then for the k-th step, we are removing the projections
of the k − 1 columns in Sk−1 selected before, now denote the singular value decomposition (V̂Sk−1 )T = USHT , then the projection matrix on columns of (V̂Sk−1 )T is UUT .
Also note that V̂(V̂Sk−1 )T = (V̂U)SHT , it is easy to check that this is singular value

33

decomposition of V̂(V̂Sk−1 )T , and the projection matrix on columns of V̂(V̂Sk−1 )T is
V̂U(V̂U)T = V̂UUT V̂T . Now the norm we need to pick from for SPA on V̂V̂T is
k(I − V̂UUT V̂T )V̂V̂T ei k = kV̂(I − UUT )V̂T ei k = k(I − UUT )V̂T ei k,
so the norms to pick for SPA on V̂T and on V̂V̂T will still be same and picking max norm
will also give the same index.
′
Lemma B.1. (Theorem 3 of Gillis et al. [12]). Let MP
= M + N = WH + N ∈ Rm×n ,
r×n
′
m×r
′
where M = WH = W[Ir |H ] ∈ R
, H ∈ R+ and rk=1 H′kj ≤ 1, ∀j and r ≥ 2. Let
K(W) = maxi kW:,i k2 , and kN:,i k2 ≤ ǫ for all i with

ǫ < σr (W) min



1
1
√
,
2 r−1 4


−1
K(W)2
1 + 80 2
σr (W)

and J be the index set of cardinality r extracted by SPA. Then there exists a permutation P
of {1, 2, · · · , r} such that


K(W)2
.
max kM′:,J(j) − W:,P (j) k ≤ ǭ = ǫ 1 + 80 2
1≤j≤r
σr (W)
Theorem B.2. Let Sp be the indices set returned by SPA in Algorithm 2, V̂p = V̂(Sp , :).
If Assumptions 2.1, 2.2, and 3.1 are satisfied, then ∃ a permutation matrix Π ∈ RK×K such
that
 √ 
V̂p − ΠT VP (VT V̂) = O ν Kǫ
F

with probability larger than 1 − O(Kn−2 ), where ǫ = Õ

error from Theorem 3.1, and rows of VP ∈ R

K×K



K 2.5 ν 1.5 (1+α0 )1.5
√
n ρλ∗ (B)



is the row-wise

are corresponding to pure nodes.

Proof of Theorem B.2. Recall we have for a ﬁxed i ∈ [n], with probability larger than
1 − Kn−3 , eTi V̂V̂T − eTi VVT = ǫ from Theorem 3.1.

Note that from Lemma 2.3, V = ΘVP . Let M′ = V̂V̂T , W = VVTP , H = ΘT , r = K,
then for M′ = WH + N, we have kN:,i k2 ≤ ǫ with probability larger than 1 − O(Kn−2 ).
W.L.O.G., let the ﬁrst K rows of Θ be K diﬀerent pure nodes. Now use Lemma B.1, there
exists a permutation π of {1, 2, · · · , K} such that


K(W)2
T
= O (νǫ) ,
max kV̂V̂p (:, j) − W(, : π(j))k = ǫ 1 + 80 2
1≤j≤K
σK (W)
√
since K(W) = maxi kW:,i k2 ≤ σ1 (W) and σσK1 (W)
(W) = κ(W) ≤ κ(VP ) = O( ν) with
probability larger than 1 − Kn−3 , from Lemma A.2.
So ∃ a permutation matrix Π ∈ RK×K such that
√
kV̂V̂Tp − WΠkF = K · O (νǫ) ,

taking transpose, it gives
 √ 
kV̂p V̂T − ΠT VP VT kF = O ν Kǫ ,

34

and
V̂p − ΠT VP (VT V̂)

F

=



V̂p V̂T − ΠT VP VT V̂

with probability larger than 1 − O(Kn−2 ).

F

 √ 
= O ν Kǫ

Theorem B.3. Let Sp be the set of of pure nodes extracted using SPACL. Let V̂p denote the
rows of V̂ indexed by Sp , and VP denote the pure nodes of V. Then, if Assumptions 2.1, 2.2,
and 3.1 are satisfied,
√
T
T
−1
kV̂−1
kF ≤ O(ν Knǫ)
p − (Π VP (V V̂))
 2.5 1.5

1.5
0)
with probability larger than 1 − O(Kn−2 ), where ǫ = Õ K nν√ρλ(1+α
is the row-wise
∗ (B)

error from Theorem 3.1, and rows of VP ∈ RK×K are corresponding to pure nodes.

Proof of Theorem B.3. In what follows we use Mp to denote the rows of M indexed by Sp
when M is n × K and by the square submatrix M (Sp , Sp ) is when M is n × n. We now
combine the row-wise deviation bound from Theorem 3.1 with Theorem B.2 to derive that,
∃ a permutation matrix Π ∈ RK×K , such that
T
T
−1
kV̂−1
kF = k(ΠT VP (VT V̂))−1 (ΠT VP (VT V̂) − V̂p )V̂−1
p − (Π VP (V V̂))
p kF

(30)

T
T
−1
≤ k(VT V̂)−1 kkV−1
P kkΠ VP (V V̂) − V̂p kF kV̂p k

th
First we will prove a bound on kV̂−1
singular value of V̂p ,
p k. Let σ̂i be the i

kV̂−1
p k =

1
.
σ̂K

(31)

p
√
From Lemma A.2, σK (VP ) = Ω(1/ n) and σ1 (VP ) = O( ν(1 + α0 )/n) with probability larger than 1 − Kn−3 .
Now for using the orthogonal matrix Ô ∈ RK×K constructed using Deﬁnition A.1,






V̂p V̂T − ΠT VP VT V̂ = V̂p − ΠT VP Ô + ΠT VP ÔV̂T − VT V̂,
by Lemma A.2, we have,

kV̂p − ΠT VP ÔkF ≤ kV̂p V̂T − ΠT VP VT kF · kV̂k + kVP k · kÔV̂T − VT kF · kV̂k
!
!
r
√
 √ 
ν(α0 + 1)
ν K(1 + α0 )
= O ν Kǫ + O
·O
√
n
ρnλ∗ (B)
 √ 
= O ν Kǫ
with probability larger than 1 − O(Kn−2 ).
(32)

Now, Weyl’s inequality for singular values gives us:

 √ 
|σ̂i − σi (VP )| ≤ kΠT VP Ô − V̂p k ≤ kΠT VP Ô − V̂p kF = O ν Kǫ
!
!!
r
r



 √
ν(1 + α0 )
νKn
1 
σ̂1 ≤ O
1+O
ǫ
.
1 − O ν Knǫ ,
σ̂K ≥ Ω √
n
1 + α0
n
(33)

35

Plugging this into Equation (31) we get:
 √

√ 
√ 
1
+
O
ν
n
Knǫ
=O n .
kV̂−1
k
=
O
p

√
The last step is true since Assumption 3.1 implies ν Knǫ = o(1). Note that we also have
kV−1
P k=
Since

√ 
1
=O n .
σK (VP )

√
8KkA − Pk
σK (V V̂) = σK (V (R + VÔ)) ≥ σK (Ô) − kV Rk ≥ 1 −
,
λ∗ (P)
T

T

(i)

T

(ii)

where (i) follows Weyl’s inequality, (ii) follows Lemma A.1.
Then we have
k(VT V̂)−1 k =

1
σK (VT V̂)

= O(1)

using Assumption 3.1 and Lemma A.2.
Finally putting everything together with Equation (30) we get:
T
T
−1
kV̂p−1 − (ΠT VP (VT V̂))−1 kF ≤ k(VT V̂)−1 kkV−1
P kkΠ VP (V V̂) − V̂p kF kV̂p k
√
√
−1
≤ O(1) · O(ν Kǫ) · kV−1
P kkV̂p k = O(ν Knǫ),

with probability larger than 1 − O(Kn−2 ).
Proof of Theorem 3.5. We break this up into proofs of Equations (1) and (2).
Proof of Equation (1). Recall that Θ̂ = V̂V̂−1
p . We have ∀i ∈ [n],




−1
eTi Θ̂ − ΘΠ = eTi V̂V̂−1
p − VVP Π

T
T
−1
≤ eTi (V̂ − V(VT V̂))V̂−1
+ eTi V(VT V̂)(V̂−1
)
p
p − (Π VP (V V̂))


√
V̂−1
+ eTi V · O(ν Knǫ)
≤ eTi V̂ − V(VT V̂)
p
r
(i)
√
√
2ν(1 + α0 )
≤ ǫ · O( n) +
· O(ν Knǫ)
n
!


p
√ 
K 2.5 ν 1.5 (1 + α0 )1.5
1.5
K(1 + α0 ) · Õ
·O n
=O ν
√ ∗
n ρλ (B)
!
2
K 3 ν 3 (1 + α0 )
= Õ
,
√
ρnλ∗ (B)

√
(i) follows from Lemma 3.2 and Theorem B.3. Also as kei Θk ≥ 1/ K, we have


!
2
eTi Θ̂ − ΘΠ
K 3.5 ν 3 (1 + α0 )
= Õ
,
√
ρnλ∗ (B)
eTi Θ

36

with probability larger than 1 − O(Kn−2 ).

Proof of Equation (2). Note ρ̂B̂ = V̂p ÊV̂Tp and ρB = VP EVTP . Note that kEk ≤
maxi keTi Pk1 = O(ρn), and kÊk ≤ kEk + kA − Pk = O(ρn) using Weyl’s inequality and
Theorem 5.2 of [19]. Then we have the following decomposition
ρ̂B̂ − ρΠT BΠ = V̂p ÊV̂Tp − ΠT VP EVTP Π
F



 F


T
T
T
= V̂p − Π VP Ô ÊV̂p + Π VP ÔÊ − EÔ V̂Tp + ΠT VP EÔ V̂Tp − ÔT VTP Π

F

V̂p + kVP k ÔÊ − EÔ
V̂p + kVP k kEk V̂Tp − ÔT VTP Π
F
F
F
!
!
!
r
r
r
 √ 
ν(1 + α0 )
ν(1 + α0 )
ν(1 + α0 )
+O
ÔÊ − EÔ O
≤2 · O ν Kǫ · O(ρn) · O
n
n
n
F




p
ν(1 + α0 )
=O ρν 3/2 (1 + α0 )Knǫ + O
ÔÊ − EÔ
n
F
≤ V̂p − ΠT VP Ô

Ê

using Equations (32) and (33) and Lemma A.2. Since kBkF = Ω(1), by Lemma A.1,





p
kρ̂B̂ − ρΠT BΠkF
ν(1 + α0 )
≤ O ν 3/2 (1 + α0 )Knǫ + O
· ÔÊ − EÔ
kρBkF
ρn
F
!


1.5


2.5
1.5
p
√ 
K ν (1 + α0 )
ν(1 + α0 )
· O K 2 nρ
= O ν 3/2 (1 + α0 )Kn · Õ
+O
√ ∗
n ρλ (B)
ρn
!
2
K 3 ν 3 (1 + α0 )
= Õ
√
∗
ρnλ (B)
with probability larger than 1 − O(Kn−2 ).

C

Proof of Lemma 5.5

Proof of Lemma 5.5. For t ≥ 2, we claim that this result follows via straightforward modiﬁcations of the proof of Lemma 7.10 in [11], where the main two elements are:
1. E[|Hij |m ] ≤ n1 for m ≥ 2. Note that for our setting, Assumption 2.1 implies that
nρ ≥ 1. Hence |Hij | ≤ 1, and hence
E[|Hij |m ] ≤ E[|Hij |2 ] ≤

Pij
1
≤
nρ
n

2. The authors use a higher order Markov inequality. This inequality upper bounds the
number of terms that are non-zero in the summand via a multigraph construction for
path counting. Then these non-zero elements are bounded by their absolute value
and hence, even though v does not equal e, just the fact that it is ﬁxed and hence
independent of Hij , is enough to apply the proof directly to get the required result.
Using an almost identical argument as [11], we have:
E[|eTi Ht v|p ] ≤ (tp)tp kvkp∞

37

Now a higher order Markov inequality, with p = (log n)ξ /2t gives:
(tp)tp kvkp∞
1
= √ (log n)ξ
(log n)tpξ kvkp∞
2
√
ξ
= exp(−(log n) log 2) ≤ exp(−(log n)ξ /3)


P |eTi Ht v| ≥ (log n)tξ kvk∞ ≤

D

Why Pruning Works

Proving the pruning algorithm requires strong distributional conditions on the residuals
of the rows of eigenvectors. Here we present a heuristic argument of why pruning works.
Note that in the pruning algorithm, essentially we are estimating the density of points in
an ǫ-ball around every point i which has suﬃciently large norm. This should work only
if the points outside the population simplex have lower density in their ǫ-balls than the
corners of the simplex. Otherwise, the pruning will remove the corners of the population
simplex, diminishing the quality of the pure nodes. We consider K ∈ {2, . . . , 10} and

Figure 5: Top panel: fraction of nodes with high norm. Bottom panel: fraction of
nodes with high norm pruned. We vary K ∈ {2, . . . , 10} on the X axis and vary
n ∈ {2000, 3000, . . . , 6000} on the Y axis.
n ∈ {2000, 3000, . . . , 6000}, α = 1K /K, Bii = 1, Bij = 0.001 and ρ = log n/n. For each
combination we use ǫ as the median of the row-wise diﬀerence of the empirical eigenvectors
from their suitably rotated population counterpart. Let y = maxi kVi k denote the largest
row-wise norm of the population eigenvectors; recall that this occurs at one of the corners of
the simplex. Let S0 denote the set of nodes with high empirical eigenvector row-norms (the
“high-norm” nodes), deﬁned as S0 := {i : kV̂i k ≥ y+ǫ}. SPA will choose at least one of these
nodes (and possibly several of them) as its estimated corners. Let B(x, ǫ) denote the ℓ2 ball
of size ǫ centered at point x. For each of the K corners ci of the population simplex (ci equals
some row of VP ), we compute the number of neighborhood points xi := |{j|V̂j ∈ B(ci , ǫ)}|;
let δ := mini xi be the minimum neighborhood size among these corners. Similarly, for each
i ∈ S0 , we compute zi = |{j|V̂j ∈ B(V̂i , ǫ)}|. Now we count the fraction of nodes in S0 that

38

could P
be pruned without pruning the corners ci of the population simplex. This fraction is
1{zi <δ}

0
. Figure 5 shows that for almost all combinations of K and n, we have
m = i∈S|S
0|
m = 1, i.e., all the nodes in S0 do get pruned, except for K = 10, n = 2000. This is expected,
since for large K and small n the pure node density around the corners of the population
simplex will be small. Figure 5 shows the fraction |S0 |/n of high-norm nodes. For all (K, n)
combinations pruning removes about a 2% to 6% of the nodes.


E

Extra simulation results

(i)

(ii)

Figure 6: (i)Relative error against increasing K. (ii) Error against ǫα : α = (0.5 −
ǫα , 0.5, 0.5 + ǫα ).
Changing K: In Figure 6 (i) we plot relative error against increasing K. We use
ρ = 0.15, αi = 1/K, Bii = 1, i ∈ [K], Bij = 0.4 for i 6= j. We see that SPACL is a little
better than GeoNMF and much better than SVI.
Changing α: In Figure 6 (ii) we use α = (0.5 − ǫα , 0.5, 0.5 + ǫα ) and plot the relative
error against ǫα . We set K = 3, ρ = 0.15, Bii = 1, i ∈ [K], Bij = 0.5 for i 6= j. Recall
that for skewed α we get unbalanced cluster sizes. SVI behaves poorly, and SPACL is much
more stable (much smaller variance) than GeoNMF.

39

