Adaptive Estimation in Structured Factor Models with
Applications to Overlapping Clustering

arXiv:1704.06977v3 [stat.ME] 22 Mar 2018

Xin Bing∗

Florentina Bunea†

Yang Ning‡

Marten Wegkamp§

Abstract
This work introduces a novel estimation method, called LOVE, of the entries and structure of a
loading matrix A in a sparse latent factor model X = AZ + E, for an observable random vector
X in Rp , with correlated unobservable factors Z ∈ RK , with K unknown, and independent
noise E. Each row of A is scaled and sparse. In order to identify the loading matrix A, we
require the existence of pure variables, which are components of X that are associated, via A,
with one and only one latent factor. Despite the fact that the number of factors K, the number
of the pure variables, and their location are all unknown, we only require a mild condition on
the covariance matrix of Z, and a minimum of only two pure variables per latent factor to show
that A is uniquely defined, up to signed permutations. Our proofs for model identifiability are
constructive, and lead to our novel estimation method of the number of factors and of the set of
pure variables, from a sample of size n of observations on X. This is the first step of our LOVE
algorithm, which is optimization-free, and has low computational complexity of order p2 . The
second step of LOVE is an easily implementable linear program that estimates A. We prove
that the resulting estimator is minimax rate optimal up to logarithmic factors in p.
The model structure is motivated by the problem of overlapping variable clustering, ubiquitous in data science. We define the population level clusters as groups of those components of
X that are associated, via the sparse matrix A, with the same unobservable latent factor, and
multi-factor association is allowed. Clusters are respectively anchored by the pure variables,
and form overlapping sub-groups of the p-dimensional random vector X. The Latent model
approach to OVErlapping clustering is reflected in the name of our algorithm, LOVE.
The third step of LOVE estimates the clusters from the support of the columns of the
estimated A. We further guarantee cluster recovery with zero false positive proportion, and
with false negative proportion control. The practical relevance of LOVE is illustrated through
the analysis of an RNA-seq data set, devoted to determining the functional annotation of genes
with unknown function.

Keywords:

Overlapping clustering, latent model, identification, high dimensional estimation, minimax
estimation, pure variables, group recovery, support recovery, sparse loading matrix, matrix factorization.
∗

Department of Statistical Science, Cornell University, Ithaca, NY. E-mail: xb43@cornell.edu.
Department of Statistical Science, Cornell University, Ithaca, NY. E-mail: fb238@cornell.edu.
‡
Department of Statistical Science, Cornell University, Ithaca, NY. E-mail: yn265@cornell.edu.
§
Department of Mathematics and Department of Statistical Science, Cornell University, Ithaca, NY. E-mail:
mhw73@cornell.edu.
†

1

1

Introduction

In this work we consider the problem of estimating the p × K sparse loading matrix A that
parametrizes the factorization of a zero-mean observable random vector, X ∈ Rp as
X = AZ + E

(1)

from n i.i.d. realizations of X. The zero mean random vector Z ∈ RK is unobservable, and can be
viewed as a latent factor vector. E ∈ Rp is a zero-mean, unobservable random noise vector, with
independent entries. The number of factors K is not known, and both p and K are allowed to grow,
and be larger, than n. Factor models have been used as dimension reduction devices in virtually
any scientific discipline for nearly a century, and generated an enormous amount of literature. We
refer to the classical monographs of Bollen [1989] and Anderson [2003] for earlier work, and to
Izenman [2008] for a more recent survey and applications.
In this work, we revisit some of the open problems in factor model definition and estimation,
and also consider one of their much less explored usage, to overlapping clustering. For the latter,
we deem two components Xi and Xj of X similar if they have non-zero association, via the matrix
A, with the same latent factor Za . Similar variables are placed in the same cluster, Ga :

Ga := j ∈ {1, . . . , p} : Aja 6= 0 ,

for each a ∈ {1, . . . , K}.

(2)

Since each Xj can be associated with multiple latent factors, the clusters will overlap. The problem
of overlapping clustering is of wide-spread interest in virtually any scientific area, for instance in
neuroscience [Craddock et al., 2012, 2013] and genetics [Jiang et al., 2004; Wiwie et al., 2015], to
give a very limited number of examples. The solutions are typically algorithmic in nature, and their
quality is assessed against a ground scientific truth or via extensive simulation studies, for instance
Krishnapuram et al. [2001]; Bezdek [2013], among many others. These problems have not received
a systematic analysis in the statistical literature and, in particular, the problem of estimating
overlapping clusters of variables, with theoretical guarantees, remains largely unexplored.
We propose model-based clustering via A. However, A cannot be uniquely defined in (1),
without further restrictions, a phenomenon well understood over six decades ago. Most notably,
Anderson and Rubin [1956] provided an in-depth analysis of this problem, and proved that in the
absence of conditions on A and C := Cov(Z), A is not identifiable in model (1). We revisit some
of these conditions here, with a view towards our application to overlapping clustering.
Using the latter as motivation, we formalize our first modeling assumption on A. We consider
models (1) in which each row of A is scaled, to avoid scale ambiguities, and sparse, to avoid the
trivial situation in which each component Xj is associated with all latent factors. Specifically, we
assume that:
PK
(i)
a=1 |Aja | ≤ 1, and each row Aj· := (Aj1 , . . . , AjK ) is sparse, for j ∈ {1, . . . , p}.
P
The inequality in (i) allows for a |Aja | = 0, which renders more flexibility to model (1), relative
P
to the more commonly used equality conditions. If a |Aja | = 0, then Xj = Ej , and Xj is not
associated with any of the latent factors, via this model. The interpretation for clustering is that
2

the corresponding Xj = Ej does not belong to any cluster given by this model, which is a desired
feature in many practical applications, including the one presented in this paper in Section 6.
Condition (i) alone cannot ensure that A in model (1) is uniquely defined, as one can still
construct an invertible matrix Q such that AZ = AQQ−1 Z, with both A and AQ satisfying (i).
Moreover, A and AQ may not have the same sparsity pattern, creating ambiguity in the cluster
definition. We introduce below two additional conditions that allow us to show, in Section 2 below,
that A is identifiable. We defer a detailed discussion of related conditions to Section 4.4.
We call the following condition the pure variable assumption. Informally, it requires the existence
of at least two pure variables Xj , which are components of X associated with one and only one
latent factor. In Section 2 we provide examples that show that if pure variables do no exist, A in
(1) is not uniquely defined.
(ii) For every a ∈ {1, . . . , K}, there exist at least two indices j ∈ {1, . . . , p} such that |Aja | = 1
and Ajb = 0 for all b 6= a.
This assumption has an immediate practical implication for variable clustering. Clusters given by
(2) are defined relative to the unobservable factor Za . A pure variable Xj is an observable proxy
of Za , and helps explain the nature of Ga . For future reference, we let I denote the index set
corresponding to pure variables. In psychology, these variables are called factorially simple items
(McDonald [1999]). A similar condition can be traced back to the econometrics literature, and an
early reference is Koopmans and Reiersol [1950], further discussed in Anderson and Rubin [1956],
who called it “zero elements in specified positions”. These works prove that (ii) corresponding to a
known set I is a sufficient condition for identifying A, for latent factors with arbitrary correlations.
However, full generality on the positive definite covariance matrix C of the latent factors comes at
the steep price of knowing I a priori, which is often unrealistic in practice. Appropriate conditions
on C that guarantee identifiability of I in (ii), in the general case when I is not known, have not
been investigated for the general model (1), to the best of our knowledge.
However, the problem of identifying A under (ii), with I unknown, has been revived more
recently, in the particular case of random vectors X and Z with only non-negative values, and when
A also has only non-negative entries. This set-up corresponds to the area of non-negative matrix
factorization (NMF), in which one studies positive matrix factorizations of the type X = AZ + E,
where the observed data X is a p × n matrix, Z is the K × n unobservable matrix of the latent
vectors values, and E is the p × n noise matrix. In this context, when E = 0, and conditioning
on Z, Donoho and Stodden [2004] was among the first to propose a condition similar to (ii), with
I unknown, coupled with appropriate conditions on Z, leading to an NMF decomposition with
unique factors. Moreover, the unique determination of I in NMF under (ii) was solved in Bittorf
et al. [2012], for known K, and for scaled NMF models, in which the columns of X, Z and A sum
up to 1, under the condition that no row of a scaled version of Z is a convex combination of the
other rows. This condition implies, in particular, that there are no duplicate (scaled) latent factors,
and similar results have been established, in the context of topic models in Arora et al. [2013].
We study here model (1) under (ii), with I and K unknown, under the restriction on the row
sums of A given in (i), but with X, Z and A allowed to have arbitrary signs. This makes the
3

direct import of conditions on Z employed in the NMF literature difficult, and we provide below
a condition on C that is appropriate for our model, and is a key ingredient in the identification
of I in Section 2 below. In this work we depart from the widely used independence assumption
between the latent factors, see for instance, Anderson and Rubin [1956]; Anderson [2003]; Bollen
[1989]; Everitt [1984]; Izenman [2008] and our discussion in Section 4.4. We allow C to be far from
a diagonal, or even a sparse matrix. Instead, we require
(iii) ∆(C) := mina6=b (Caa ∧ Cbb − |Cab |) > 0,
where a ∧ b := min(a, b). If (iii) holds, then Cov(Za ± Zb ) = Var(Za ) + Var(Zb ) ± 2 · Cov(Za , Zb ) ≥
Caa + Cbb − 2|Cab | > 0, which implies that the latent factors are different, up to signs, that is
|Za | =
6 |Zb | a.s. for any a 6= b. Thus, (iii) can be viewed as a minimal natural requirement for the
study of this problem.
We note that under the pure variable assumption (ii), if I is known in advance, the arguments
employed in the proof of Theorem 2 of Section 2 show that (iii) is not required, consistent with
the existing literature on general latent models. When I is not known, our condition (iii) implies
that we do not have duplicate latent factors, up to signs, consistent with similar requirements in
NMF, as explained above. Therefore, condition (iii) can be viewed as bridging the gap between the
classical literature in general latent models and the more recent literature devoted to the specialized
NMF-type models.
Summarizing, this work is devoted to estimation in model (1) with A satisfying (i) - (ii), and
with the positive definite matrix C satisfying (iii). The number of factors K is not known, and both
K and p are allowed to grow and be larger than n. In Section 1.1 below we present our contributions
and the structure of this paper. A detailed contrast with existing literature is presented in Section
4.4.

1.1

Our contributions

1. Identifiability of the allocation matrix A in sparse latent models with pure variables.
We show, in Proposition 2 of Section 2, that the allocation matrix A, which is allowed to have entries
of arbitrary signs, is uniquely defined, up to trivial orthogonal transformations, namely signed
permutation matrices. This is a consequence of one of our main results, Theorem 1 of Section 2.
In this result we highlight and resolve the main difficulty in this problem, that of distinguishing
between the pure variables and the non-pure variables. Both proofs are constructive, and show that
the pure variable set I and allocation matrix A can be determined uniquely from Σ := Cov(X).
Moreover, the number of factors K is not assumed to be known, and its determination is also a
consequence of Theorem 1. To the best of our knowledge, these are new results in both the latent
factors literature and other related matrix factorization literature. We comment on connections to
related results in Section 4.4.
2. Estimation of the allocation matrix A and of the overlapping clusters. The LOVE
b of the sparse and structured matrix A that is tailored to our
algorithm. We provide an estimator A
model specifications. Our approach follows the constructive techniques used in our identifiability
b an estimator of the pure variable set I, and K,
b an estimator of the
proofs. We first construct I,
4

number of clusters, K. These are used to estimate the rows in A corresponding to pure variables.
The remaining rows of A are estimated via an easily implementable linear program that is tailored
to this problem. As part of our procedure, we also develop a novel estimator (12) and (13) of a
precision matrix, C −1 . Our procedure is presented in Sections 3.1, 3.2 and 3.3, respectively. To the
best of our knowledge, our estimation strategy is new, and complements the large body of literature
in factor models. In particular, we do not resort to optimizing a complicated quasi likelihood
function via computationally demanding EM algorithms. These algorithms require, in addition, a
notoriously delicate initialization, especially in high dimensions, and typically only convergence to
a stationary point can be guaranteed, see Rubin and Thayer [1982]. Moreover, as our procedure is
not Bayesian, we do not employ distributional assumptions to construct our estimator. In Section
b using the estimated allocation matrix A.
b The
3.4, we build a collection of overlapping clusters G,
combined procedure is summarized in a new algorithm, LOVE, highlighting our Latent model
approach to OVErlapping clustering.
3. Statistical guarantees. Our estimation procedure does not depend on distributional assumptions, but for the purpose of our statistical analysis, and in particular our minimax analysis, we
assume that X has a sub-Gaussian distribution. LOVE, for appropriate choices of tuning parameters, recovers the population level clusters with a zero false positive proportion and generally low
false negative proportion, with high probability, and under a mild condition on the cluster separation as measured by the quantity ∆(C). This is a direct consequence of a number of results
regarding estimation of identifiable loading matrices in factor models satisfying (i) - (iii) and, to
the best of our knowledge, are all new.
(1) Consistent estimation of the number of factors K;
(2) Control of the relationship between Ib and I for A with entries of arbitrary strength. In
particular, we show I ⊆ Ib ⊆ I ∪ J1 , where we carefully define and characterize J1 as the set
of quasi-pure variables.
b A), defined below, for all q ≥ 1, in particular for
(3) Minimax lower bounds on the norms Lq (A,
q = +∞, for A given by model (1) under (i) - (iii).

(4) Attainment of these bounds, showing that our procedure is minimax optimal and adaptive.
b
(5) Control of the relationship between the support of A and the support of A.

(6) Control of cluster recovery.

The details are given in Sections 4.2 and 4.3. In particular, we emphasize that (2) above, proved
in Theorem 3 of Section 4.1, guarantees recovery of I with minimal mistakes. This result does not
require the necessary, yet unpleasant, signal strength restrictions encountered in the typical exact
support recovery literature. However, under such restrictions, we also obtain Ib = I, with high
probability, in Remark 2 of Section 4.1. Since placing restrictions on the entries in A reduces the
number of configurations of interest, the more general result (2) is a new and practically relevant
result for pure variable recovery.
Results (3) and (4) are given in Theorems 4 and 5, respectively, of Section 4.2. We consider the

5

loss function
b − Ak∞,q ,
b A) := min kAP
Lq (A,
P

1 ≤ q ≤ ∞,

with the minimum taken over all K × K signed permutation matrices P and


kAk∞,q := max kAi· kq = max 
1≤i≤p

1≤i≤p

K
X
j=1

1/q

|Aij |q 

,

is the maximum `q norm of the rows of A. We let s = maxi∈[p] kAi· k0 be the row-sparsity index.
As explained in Section 3.3, our estimation procedure is closely connected to sparse estimation
in errors-in-variable regression models with design matrix C, and the error rates we obtain can
be interpreted in that context. We show that the error of estimation with respect to the Lq loss
function, for each q, is proportional to s1/q n−1/2 , multiplied by kC −1 k∞,1 . This is consistent with
the most recent results regarding error rates expressed in terms of the `q -sensitivity of C in Belloni
et al. [2017], as discussed in Section 4.2. The results hold up to logarithmic factors in p and s.
Results (5) and (6) are presented in Theorem 6 of Section 4.3. Moreover, we can further partition
the variables in each cluster into two signed sub-groups consistently. In our model formulation,
A is allowed to have positive and negative entries. Since A can only be identified up to signed
b However, we can identify consistently
permutations, one cannot expect sign consistency for A.
the two sub-groups of each cluster that contain variables that are associated with the common
latent factor in the same direction, although the direction itself is not identifiable. These results
are presented in Section 4.3.
We conduct an extensive simulation study in Section 5 to assess the numerical performance of
our proposed strategy. The study confirms our theoretical findings. We conclude the validation
of our approach with a data analysis, devoted to determining the functional annotation of genes
with unknown function. Our analysis confirms existing biological ground truths, as our procedure
tends to cluster together genes with the same Gene Ontology (GO) biological process, molecular
function, or cellular component terms.
We summarize our contributions in the table below, restricting attention to estimation in general
latent models (1) under (i) - (iii), without any further restrictions on the signs or scales of X and
Z.
Model (1) under (i) - (iii)
Identifiability Conditions
Estimation: I
Estimation: A
Guarantees: I
Guarantees: A
Cluster recovery

Our results
Existence of I; I and K unknown; mild conditions on C.
Optimization-free. Runs in O(p2 ) time
Not MLE-based. Uses linear program. Runs in O(p2 + pK).
Recovered
Finite sample lower bounds for k k∞,q recovery
Adaptive finite sample upper bounds
p and K can grow with n
Guaranteed

Existing results in factor models
Existence of known I; general C
×
E-M type, with known I. Computationally involved. Not guaranteed to estimate A
×
×
Row-wise asymptotic normality of MLE-type estimators
p can grow with n; K fixed.
×

In Section 4.4 we discuss our results further, and provide a detailed comparison between our
work and related contributions. All proofs are deferred to Section A of the supplementary material.

6

1.2

Notation

We use the following notation throughout this paper. For the n consecutive integer set starting
from 1, we write [n] = {1, . . . , n}. The sign of any generic number N is denoted by sign(N ). For
any m × d matrix M and index sets I ⊆ {1, . . . , m} and J ⊆ {1, . . . , d}, we write MI to denote
the |I| × d submatrix (Mij )i∈I,1≤j≤d of M consisting of the rows in the index set I, while we denote by MIJ the |I| × |J| submatrix with entries Mij , i ∈ I and j ∈ J. The ith row of M is
denoted by Mi· , and the jth column of M is denoted by M·j . Let kM k∞ = max1≤j≤m,1≤k≤d |Mjk |,
Pd
Pd
P
P
2 )1/2 , kM k
Mjk
kM k1 = 1≤j≤m,1≤k≤d |Mjk |, kM kF = ( m
∞,1 = max1≤j≤m
k=1 |Mjk |
k=1
j=1
Pm
and kM k1,∞ = max1≤k≤d j=1 |Mjk | denote the matrix max norm, matrix `1 norm, matrix Frobenius norm, matrix 1 norm and matrix ∞ norm. We denote by h · i the Frobenius scalar product.
P
For a vector v ∈ Rd , define kvkq = ( di=1 |vj |q )1/q for 1 ≤ q < ∞, kvk∞ = max1≤j≤d |vj | and
kvk0 = |supp(v)|, where supp(v) = {j : vj 6= 0} and |A| is the cardinality of the set A. We write
M T for the transpose of M and diag(m1 , . . . , md ) for the d × d diagonal matrix with elements
m1 , . . . , md on its diagonal, while diag(M ) is the diagonal matrix obtained from the diagonal elements of a square matrix M . We use c0 , c1 , . . . to denote generic constants. Finally, a signed
permutation matrix is an orthogonal matrix that permutes the index and switches the sign within
each column. We write HK as the hyperoctahedral group of K × K signed permutation matrices.

2

Identifiability

In this section we show that the allocation matrix A given by Model (1) and (i) - (iii) is identifiable,
up to multiplication with a signed permutation matrix.
For any A ∈ Rp×K which satisfies Model (1), we can partition the set [p] = {1, . . . , p} into two
disjoint parts: I and its complement J := [p] \ I such that for each row Ai· of AI , there exists only
one a ∈ [K] such that |Aia | = 1. We name I the pure variable set and J the non-pure variable set.
Specifically, for any given A, the pure variable set I is defined as
[
I(A) :=
Ia ,
Ia := {i ∈ [p] : |Aia | = 1, Aib = 0, for any b 6= a} .
(3)
a∈[K]

We write I(A) in (3) to emphasize that the pure variable set is defined relative to A. In the
following, we will not write this explicitly when there is no confusion. We also note that the sets
{Ia }1≤a≤K form a partition of I.
To show the identifiability of A, it suffices to show that AI and AJ are identifiable, respectively,
up to signed permutation matrices. By the definition of AI , this matrix is identifiable provided the
partition of the pure variable set I is. The identifiability of I, and thus the problem of distinguishing
between the sets I and J, on the basis of the distribution of X alone, is the central challenge in
this problem. We meet this challenge in Theorem 1 below: part (a) offers a necessary and sufficient
characterization of I; part (b) shows that, as a consequence, I and its partition I := {Ia }1≤a≤K
are identifiable. Let
Mi := max |Σij |
(4)
j∈[p]\{i}

7

be the largest absolute value of the entries of row i of Σ excluding |Σii |. Let Si be the set of indices
for which Mi is attained:

Si := j ∈ [p] \ {i} : |Σij | = Mi .
(5)
Theorem 1. Assume that model (1) and (i) - (iii) hold. Then:
(a) i ∈ I

⇐⇒

Mi = Mj for all j ∈ Si .

(b) The pure variable set I can be determined uniquely from Σ := Cov(X). Moreover, its partition
I := {Ia }1≤a≤K is unique and can be determined from Σ up to label permutations.
The identifiability of the allocation matrix A and that of the collection of clusters G = {G1 , . . . , Gk }
in (2) use the results from Theorem 1 in crucial ways. We state the result in Theorem 2 below.
Theorem 2. Assume that Model (1) with (i) - (iii) holds. Then, there exists a unique matrix A,
up to a signed permutation, such that X = AZ + E. This implies that the associated overlapping
clusters Ga , for 1 ≤ a ≤ K, are identifiable, up to label switching.
Remark 1. We show below that the pure variable assumption (ii) is needed for the identifiability
of A, up to a signed permutation. Assume that X = AZ + E satisfies (i) and (iii), but not (ii). We
construct an example in which X can also be written as X = ÃZ̃ + E, where Ã and Z̃ satisfy the
same conditions (i) and (iii), respectively, but Ã 6= AP for any K × K signed permutation matrix
P and Ã may have a sparsity pattern different from A. To this end, we construct Ã and Z̃ such
that ÃZ̃ = AZ. Let Ã = AQ and Z̃ = Q−1 Z, for some K × K invertible matrix Q to be chosen. By
the definition of C, Q needs to satisfy C = QCQT . In addition, we need to guarantee that Ã = AQ
satisfies (i). For simplicity, we set K = 3. The following example satisfies all our requirements:




p
1 0 0
1/2
3/8 0


 p

C = 0 2 0 , Q = −2 3/8 1/2 0 .
0 0 3
0
0
1
It is easy to verify that C = QCQT holds. For any 1 ≤ j ≤ p, consider
r
r
1 1 2 3 1 2
Aj· = ( ,
, −
)
4 8 3 4 8 3
then
1
Ãj· = Aj· Q = (0,
8

r

3
1
+
2 16

r

2 3 1
, −
3 4 8

r

2
),
3

which also satisfies condition (i). However, Aj· and Ãj· have different sparsity patterns. Thus, if
the matrix A does not satisfy (ii), A is generally not identifiable.

3

Estimation

We develop estimators from the observed data, which is assumed to be a sample of n i.i.d. copies
X (1) , . . . , X (n) of X ∈ Rp , where p is allowed to be larger than n. Our estimation procedure consists
8

of the following four steps:
(1) Estimate the pure variable set I, the number of clusters K and the partition I;
(2) Estimate AI , the submatrix of A with rows Ai· that correspond to i ∈ I;
(3) Estimate AJ , the submatrix of A with rows Aj· that correspond to j ∈ J;
(4) Estimate the overlapping clusters G = {G1 , . . . , GK }.

3.1

Estimation of I and I

Given the different nature of their entries, we estimate the submatrices AI and AJ separately. For
the former, we first estimate I and its partition I = {I1 , . . . , IK }, which can be both uniquely
constructed from Σ, as shown by Theorem 1. We use the constructive proof of Theorem 1 for this
step, replacing the unknown Σ by the sample covariance matrix
n

X
b= 1
Σ
X (i) (X (i) )T .
n
i=1

Specifically, we iterate through the index set {1, 2, . . . , p}, and use the sample version of part (a)
of Theorem 1 to decide whether an index i is pure. If it is not deemed to be pure, we add it to the
set that estimates J. Otherwise, we retain the estimated index set Sbi of Si defined in (5), which
corresponds to an estimator of Mi given by (4). We then use the constructive proof of part (b) of
Theorem 1 to declare Sbi ∪ {i} := Ib(i) as an estimator of one of the partition sets of I. The resulting
procedure has complexity O(p2 ), and we give all the specifics in Algorithm 1 of Section 3.5. The
algorithm requires the specification of a tuning parameter δ, which will be discussed in Section 5.1.

3.2

Estimation of the allocation submatrix AI

b K
b and Ib = {Ib1 , . . . , Ibb } from Algorithm 1, we estimate the matrix AI by
Given the estimators I,
K
b ×K
b matrix with rows i ∈ Ib consisting of K
b − 1 zeros and one entry equal to either +1 or −1
a |I|
b
as follows. For each a ∈ [K],

bia = 1. Note that A
bia can only be +1 or −1 by
(1) Pick an element i ∈ Iba at random, and set A
the definition of a pure variable.
bja = sign(Σ
b ij ).
(2) For the remaining j ∈ Iba \ {i}, we set A

This procedure induces a partition of Iba = Iba1 ∪ Iba2 , where Iba1 and Iba2 are defined below:
(
bka = A
bla , for k, l ∈ Iba1 or k, l ∈ Iba2
A
.
bka 6= A
bla , for k ∈ Ib1 and l ∈ Ib2
A
a

3.3

(6)

a

Estimation of the allocation submatrix AJ

We continue by estimating the matrix AJ , row by row. To motivate our procedure, we begin by
highlighting the structure of each row Aj· of AJ , for j ∈ J. We recall that Aj· is sparse, with
kAj· k1 ≤ 1, for each j ∈ J, as specified by assumption (i). In addition, model (1) subsumes a
9

further constraint on each row Aj· of A, as explained below. To facilitate notation, we rearrange Σ
and A as follows:
"
#
" #
ΣII ΣIJ
AI
Σ=
and A =
.
ΣJI ΣJJ
AJ
Model (1) implies the following decomposition of the covariance matrix of X:
"
# "
# "
#
ΣII ΣIJ
AI CATI AI CATJ
ΓII
0
Σ=
=
+
.
ΣJI ΣJJ
AJ CATI AJ CATJ
0 ΓJJ
In particular, ΣIJ = AI CATJ . Thus, for each i ∈ Ia with some a ∈ [K] and j ∈ J, we have
Aia Σij =

A2ia

K
X

Ajb Cab =

b=1

K
X

T
Ajb Cab = Ca·
Aj· .

(7)

b=1

Averaging display (7) over all i ∈ Ia yields

1 X
T
Aia Σij = Ca·
Aj· ,
|Ia |
i∈Ia

for each a ∈ [K].

(8)

For each j ∈ J, we let β j = Aj· and


T
X
X
1
1
θj = 
AiK Σij  .
Ai1 Σij , . . . ,
|I1 |
|IK |
i∈I1

(9)

i∈IK

Since Aia ∈ {−1, 1}, for each i ∈ Ia and a ∈ [K], the entries of θj are respective averages of the
sign corrected entries of Σ corresponding to the partition of the pure variable set. Summarizing,
modeling assumption (i) and equation (8) above show that the estimation of AJ reduces to estimating, for each j ∈ J, a K-dimensional vector β j that is sparse, with norm kβ j k1 ≤ 1, and that
satisfies the equation θj = Cβ j . Both C and θj , for each j ∈ J, can be estimated directly from the
b we estimate the a-th entry of θj by
data as follows. For each j ∈ J,
1 X b b
b
Aia Σij , a ∈ [K],
θbaj =
|Iba |

(10)

i∈Iba

and compute
baa =
C

1
|Iba |(|Iba | − 1)

X

i,j∈Iba ,i6=j

b ij |,
|Σ

bab =
C

1
|Iba ||Ibb |

X

i∈Iba ,j∈Ibb

bia A
bib Σ
b ij ,
A

(11)

b and b ∈ [K]
b \ {a} to form the estimator C
b of C. The estimates (10) and (11) rely
for each a ∈ [K]
crucially on having first estimated the pure variables and their partition, according to the steps
described in Sections 3.1 and 3.2 above.
We could cast the estimation of β j as a sparse, errors in variables regression problem by viewing
b as observed quantities in the model given by θbj = Cβ j +ε and C
b = C +η, with “unobserved
θbj and C
10

b − C. Although the “errors” ε and η are heavily correlated in this
errors” ε = θbj − θj and η = C
context, unlike the typical assumptions on the errors in such models, it is nevertheless tempting
to use an appropriate extension of the optimal estimation procedure developed for sparse errors in
variables regression in Belloni et al. [2017]. This would essentially boil down to minimizing kβ j k1
b j − θbj k∞ . While Belloni et al. [2017] showed that their conic
under a suitable constraint on kCβ
programming estimator is optimal in a minimax sense, conic programming for each row j ∈ Jb would
be computationally costly and involve too many tuning parameters. Moreover, our final goal of
support recovery of β j would require additional thresholding for each row. The latter drawback is
explicitly mentioned in the discussion following Theorem 1 in Belloni et al. [2017].
For these reasons, we do not pursue this line of estimation here. Instead, we have developed a
more computationally efficient method, that only involves two tuning parameters, but shares the
same rate of convergence of the resulting estimators as the aforementioned procedure in Belloni
et al. [2017]. We exploit the fact that the square matrix C is invertible and take the equation
b θbj , based
β j = C −1 θj as our starting point. The idea is to first construct a pre-estimator β̄ j = Ω
−1
b of the precision matrix Ω := C , followed by a sparse projection
on an appropriate estimator Ω
j
of β̄ . Alternatively, we could use a simple hard threshold in the second step.
We first motivate our proposed estimator of Ω. From the decomposition
β̄ j − β j
we immediately have

b θbj − θj ) + (Ω
b − Ω)θj
= Ω(
b θbj − θj ) + (ΩC
b − I)β j ,
= Ω(

b ∞,1 kθbj − θj k∞ + kΩC
b − Ik∞ kβ j k1 .
kβ̄ j − β j k∞ ≤ kΩk

(12)

(13)

Since we can show, in Lemma 11, that kθbj −θj k∞ has optimal convergence rate, and since kβ j k1 ≤ 1
b should ideally render values for kΩk
b ∞,1 and kΩC
b − Ik∞ that are
under our model, our estimator Ω
as small as possible. With this in mind, we propose the linear program

subject to

b t̂) = arg
(Ω,

Ω = ΩT ,

min

b K
b
t∈R+ , Ω∈RK×

b − Ik∞ ≤ λt,
kΩC

t

kΩk∞,1 ≤ t,

(14)

(15)

with tuning parameter λ. This linear programming problem clearly is tailored to our purpose, and
b adds a novel estimator for C −1 to the rich literature of precision matrix
its optimal solution Ω
estimation (Meinshausen and Bühlmann [2006]; Yuan and Lin [2007]; Friedman et al. [2008]; Cai
et al. [2011, 2016], to name a few). Its novelty consists in (a) the usage of the matrix k · k∞,1 norm,
rather than the commonly used matrix k · k1 norm, and (b) the fact that this norm appears in the
upper bound of the restriction (15). Our analysis does not require any additional assumptions on
C, i.e., no sparsity assumption on C or C −1 is imposed. In particular, we avoid the common, albeit
highly undesirable, assumption that the condition number of C is bounded [Cai et al., 2016]. After
b θbj , for each j ∈ J,
b we solve the following optimization problem
we compute β̄ j = Ω
βbj = arg min kβk1
b
β∈RK

11

(16)

subject to
kβ − β̄ j k∞ ≤ µ,

(17)

for some tuning parameter µ that is constant in j, to obtain our final estimate βbj as the optimal
solution of this linear program. This solution is sparse and properly scaled, in accordance to our
b we obtain A
b b, and merging A
bb with
model specification (i). By stacking βbj over all rows j ∈ J,
J
I
b b produces our final estimator A
b of A. Its statistical property is analyzed in Section 4 along with
A
J
precise forms of the tuning parameters.

3.4

Estimation of the overlapping groups

Recalling the definition of groups in (2), the overlapping groups are estimated by


b1 , . . . , G
bb , G
b a = i ∈ [p] : A
bia 6= 0 , for each a ∈ [K].
b
Gb = G
K

(18)

as the pure noise cluster. We can estimate G0 in (19) by

b 0 = i ∈ [p] : A
bia = 0, for all a ∈ [K]
b .
G

(20)

b with the same latent factor Za are therefore placed in the
Variables Xi that are associated (via A)
b a . To accommodate potential pure noise variables, we further define
same group G

G0 := j ∈ {1, . . . , p} : Aja = 0, for all a ∈ {1, . . . , K}
(19)

However, our main focus is on G because it completely determines G0 .
In many applications, it may be of interest to identify the sub-groups of variables that are all
either positively or negatively associated with the same latent factor. To this end, we define
n
o



G s := Gs1 , . . . , GsK , Gsa := G1a , G2a :=
(21)
i ∈ Ga : Aia > 0 , i ∈ Ga : Aia < 0 ,

and they are estimated by
 s
b1 , . . . , G
bs ,
Gbs = G
b
K
n
o

s
ba =
b
ba : A
bia > 0 , i ∈ G
ba : A
bia < 0 , for each a ∈ [K].
G
i∈G

(22)

The fact that A is only identifiable up to a signed permutation matrix, has the repercussion that
the labels of the two sub-groups in Gsa are not identifiable. Thus, variables placed in the subgroups
G1a and G2a are, respectively, associated with Za in the same direction. The directions between two
sub-groups, henceforth called direction sub-groups, are opposite. This can be identified, although
the direction itself cannot. We show in Section 4 that the direction sub-groups can be identified,
and well estimated.

3.5

LOVE: A Latent variable model approach for OVErlapping clustering.

We give below the specifics of Algorithm 1, motivated in Section 3.1, and summarize our final
algorithm, LOVE in Algorithm 2.
12

Algorithm 1 Estimate the partition of the pure variables I by Ib
b δ)
1: procedure PureVar(Σ,
b
2:
I ← ∅.
3:
for all i ∈ [p] do

b ij | ≤ |Σ
b il | + 2δ
4:
Ib(i) ← l ∈ [p] \ {i} : maxj∈[p]\{i} |Σ
5:
P ure(i) ← T rue.
6:
for all j ∈ Ib(i) do
b ij | − maxk∈[p]\{j} |Σ
b jk | > 2δ then
7:
if |Σ
8:
P ure(i) ← F alse,
9:
break
10:
end if
11:
end for
12:
if P ure(i) then
13:
Ib(i) ← Ib(i) ∪ {i}
b
14:
Ib ← Merge(Ib(i) , I)
15:
end if
16:
end for
b
17:
return I.
18: end procedure

b
19: function Merge(Ib(i) , I)
20:
for all G ∈ Ib do
21:
if G ∩ Ib(i) 6= ∅ then
22:
G ← G ∩ Ib(i)
23:
return Ib
24:
end if
25:
end for
26:
Ib(i) ∈ Ib
27:
return Ib
28: end function

4

. Ib is a collection of sets

. Replace G ∈ Ib by G ∩ Ib(i)

. add Ib(i) in Ib

Statistical guarantees

We provide in this section statistical guarantees for:
b
(1a) The estimated number of clusters K;

b
(1b) The estimated pure variable set Ib and its estimated partition I;
(2)

(3)

b and its adaptation to the unknown row sparsity of A.
The estimated allocation matrix A

The individual Group False Positive Proportion (GFPP), the individual Group False Negative
Proportion(GFNP), the Total False Positive Proportion (TFPP) and the Total False Negative
Proportion (TFNP) for the estimated overlapping groups.

13

Algorithm 2 The LOVE procedure for overlapping clustering.
b from I.I.D. data (X (1) , ..., X (n) ), the tuning parameters δ, λ and µ.
Require: Σ
b the estimated set of pure variables Ib and its
1: Apply Algorithm 1 to obtain the number of clusters K,
b
partition of I.
bb from (6).
2: Estimate AI by A
I
−1
b from (14) and β̄ j for each j ∈ J.
b
3: Estimate C
by Ω
b
b
b
b
4: Estimate AJ by AJb from (16). Combine AIb with AJb to obtain A.
b 1 , ..., G
b } and its direction subgroups Gbs = {G
b s , ..., G
b s } from (18)
5: Estimate overlapping groups Gb = {G
1

K̂

b
- (22) by using A.
b Gb and Gbs .
6: Output A,

K̂

We make the blanket assumption for the remainder of this paper that X is sub-Gaussian, that
is, the Orlicz norm kXj kψ2 of each Xj is bounded by a common constant σ.1 The sub-Gaussian
condition implies maxj∈[p] Σjj ≤ 2σ and kCk∞ ≤ 2σ. Let
E = E(δ) :=



max

1≤i<j≤p


b
Σij − Σij ≤ δ .

(23)

p
We assume throughout that δ = c0 log(p ∨ n)/n, for some constant c0 , and log p = o(n), so that
δ = o(1), for n large enough, where a ∨ b = max(a, b). Taking c0 = c0 (σ) > 0 large enough, Lemma
2 in Bien et al. [2016] guarantees that E holds with high probability:
P(E) ≥ 1 − c1 (p ∨ n)−c2

(24)

for some positive, finite constants c1 and c2 . Apart from δ, the quantity
∆(C) := ν > 0,

(25)

plays an important role in our analysis. Indeed, assumption (iii) requires that ν > 0 in order to
guarantee that the latent factors are distinguishable from one another. We can view ν as a measure
of their separation, and naturally therefore, the size of ν impacts the quality of all our estimators,
in addition to the magnitude of δ.

4.1

b Ib and Ib
Statistical guarantees for K,

We first analyze the performance of our estimator Ib of I, and its corresponding partition. This
problem belongs to the general class of pattern recovery problems, and it is well understood that
under strong enough signal conditions one can expect Ib = I, with high probability. This turns out
to be indeed the case for our problem, but we obtain this as a corollary of a more general result. We
set out to quantify when our estimated set contains the least taxing type of errors, under minimal
assumptions. To make this precise, we introduce the concept of quasi-pure variables. A quasi-pure
variable Xi has very strong association with only one latent factor, say Za , in that |Aia | ≈ 1, and
1

The Orlicz norm of Xj is defined as kXj kψ2 = inf {c > 0 : E [ψ2 (|Xj |/c)] < 1} , based on the Young function
ψ2 (x) = exp(x2 ) − 1.

14

very low association with the rest: |Aib | ≈ 0, for all b 6= a. Formally, we define the set of quasi-pure
variables as:
J1 := {j ∈ J : there exists a ∈ [K], such that |Aja | ≥ 1 − 4δ/ν}.

(26)

For each a ∈ [K] we further define the set of quasi-pure variables associated with the same factor:
J1a := {j ∈ J1 : |Aja | ≥ 1 − 4δ/ν}.

(27)

When ν is a strictly positive constant,  := 4δ/ν = o(1). The lower bound |Aja | ≥ 1 −  in (27)
implies, under condition (ii), that |Ajb | ≤ , for any b 6= a and j ∈ J1a , justifying the name quasipure variables for those components of X with indices in J1 . We observe, for future reference, that
{J11 , . . . , J1K } forms a partition of J1 .
We show in Theorem 3 that, with very high probability, the estimated Ib contains the pure
variable set I, and is in turn contained in a set that includes all pure variables and quasi-pure
variables. Importantly, Ib will not include indices of variables Xj that are associated with multiple
latent factors at a level higher than . Equally importantly, if a quasi-pure variable Xi is included
b then this variable will have the corresponding |Aia | ≈ 1, and it will be placed together with
in I,
the pure variables associated with the same factor Za , for some a, and not in a new cluster. This is
crucial for ensuring that the number of clusters K is consistently estimated, and also for establishing
the cluster misclassification proportion in Section 4.3 below.
Theorem 3. Assume Model (1) with (i) - (iii), and

 p
ν > 2 max 2δ, 2kCk∞ δ .

(28)

Then:

b = K;
(a) K

(b) I ⊆ Ib ⊆ I ∪ J1 .

Moreover, there exists a label permutation π of the set {1, . . . , K}, such that the output Ib =

Iba a∈[K] from Algorithm 1 satisfies:
π(a)
(c) Iπ(a) ⊆ Iba ⊆ Iπ(a) ∪ J1 .

All results hold with probability larger than 1 − c1 (n ∨ p)−c2 , for c1 , c2 positive constants defined in
(24).
The conclusion of Theorem 3 holds only under condition (28), which stipulates that the separation
between the latent factors, as measured by ν, is not only strictly positive, which was needed for
identifiability, but slightly above a quantity that depends on the estimation error δ, and which
becomes o(1) for n large enough.

15

Remark 2. Let e1 = (1, 0, . . . , 0)T and HK be the hyperoctahedral group of signed permutation
matrices. If AI and AJ are well separated in the sense that
min

j∈J, P ∈HK

Aj· − P e1

1

> 4δ/ν,

then J1 = ∅, and Theorem 3 yields exact recovery of the pure variable set and of its partition:
Ib = I and Ib = I, with high probability. However, we expect J1 6= ∅, as we expect quasi-pure
variables to be present in a high dimensional model, which is the context for which Theorem 3 has
been established.

4.2

b
Statistical guarantee for A

b obtained in
In this section we state, and comment on, the statistical properties of the estimate A
p
Sections 3.2 and 3.3. Recall that δ = O( log(p ∨ n)/n) was given in (23) above, and the estimation
of AJ made use of two tuning parameters: λ, in (15), and µ, in (17). Theorem 4 establishes the
properties of our estimates relative to the theoretically optimal values of these tuning parameters,
both of which are functions of δ, while their data adaptive calibration is discussed in Section 5.1
below. We let λ = 2δ 0 and µ = 5kC −1 k∞,1 δ 0 , with


8
δ0 =
kCk∞ − 3 δ,
(29)
ν
for ν defined in (25) above. When ν and kCk∞ are strictly positive constants we thus have
p
p
λ = O( log(p ∨ n)/n) and µ = O(kC −1 k∞,1 log(p ∨ n)/n). We consider the loss function for two
p × K matrices A, A0 as
Lq (A, A0 ) := min kAP − A0 k∞,q ,

1 ≤ q ≤ ∞.

P ∈HK

(30)

Here HK is the hyperoctahedral group of all K × K signed permutation matrices and


kAk∞,q := max kAi· kq = max 
1≤i≤p

1≤i≤p

K
X
j=1

for a generic matrix A ∈ Rp×K .

1/q

|Aij |q 

,

Theorem 4. Assume the conditions in Theorem 3 hold. Let λ and µ be as defined above, and set
s = maxi∈[p] kAi· k0 . Then,
b A) ≤ 10s1/q kC −1 k∞,1 δ 0 ,
Lq (A,

1 ≤ q ≤ ∞,

with probability larger than 1 − c1 (n ∨ p)−c2 , for c1 , c2 positive constants defined in (24), provided
that (2µ + 4δ/ν) < 1. We use the convention that s1/q = 1 for q = +∞.
Remark 3.

16

1. In fact, we prove the stronger result
bi· − (AP )i·
min A

P ∈HK

q

≤ 10(si )1/q kC −1 k∞,1 δ 0 ,

1 ≤ q ≤ ∞,

with sparsity index si = kAi· k0 for each row Ai· , i ∈ [p] of A. The signed permutation matrix
P that achieves the minimum is determined by the alignment of the pure variables and is the
same for each i ∈ [p].
bi· k1 ≤ 1, for each i ∈ [p], with high
2. Inspection of the proof of this result quickly reveals that kA
probability, in accordance with our model requirement (i).
3. Our primary focus is the bound for q = +∞, as this leads to inference on support recovery
of A. More generally, for any q ≥ 1, it is well understood that the quality of estimating a
sparse vector in high-dimensional regression-type models depends on the interplay between its
sparsity and the behavior of the appropriate Gram matrix associated with the model, which
reduces to C = E[ZZ T ] in our case. The concept of `q -sensitivity, introduced by Gautier and
Tsybakov [2011]; Belloni et al. [2017], is the most general characterization of this interplay
to date. It facilitates a link between the `q -norm of sparse vectors β and the `∞ -norm of the
product between the Gram matrix and β, uniformly over vectors β of sparsity s, ranging over
a collection of cones. In our context, that of a square, invertible matrix C, the reciprocal of
the `∞ -sensitivity of C becomes essentially kC −1 k∞,1 , which indeed links kβk∞ to kCβk∞ .
Similarly, the quantities s1/q kC −1 k∞,1 provide concrete substitutes of the reciprocals of the
`q -sensitivities of C, and all of our rates in Theorem 4 match the lower bounds in Theorem
5, up to a logarithmic factor in p/s and the quantity kC −1 k∞,1 .

The size of kC −1 k∞,1 ranges from the constant kC −1 k∞ , when all latent factors are independent, to the fully general case of kC −1 k∞,1 = O(K). In the latter case the bounds become
p
meaningful when K < O( n/ log p). However, if C −1 is sparse, then kC −1 k∞,1 may be considerably smaller than K. In particular, if Z has a multivariate normal distribution and many
factors Zi are conditionally independent, then kC −1 k∞,1 is small. We do not make any of
these assumptions here, and regardless of the situation, Theorem 4 shows that our estimation
procedure adapts automatically to it.

We now show that the rates of convergence in Theorem 4 are optimal (up to a logarithmic
factor in p/s and provided kC −1 k∞,1 remains bounded) in a minimax sense for all estimators over
the parameter space


As := A ∈ [−1, 1]p×K : A satisfies conditions (i) and (ii) and max kAi· k0 ≤ s .
1≤i≤p

For our purpose of establishing a minimax lower bound, it suffices to consider a particular subGaussian distribution of (Z, E). We choose to take the multivariate Gaussian with E[(Z, E)] = 0 ∈
RK+p , Cov(Z, E) = 0 ∈ RK×p and Cov(Z) = IK , the K × K identity matrix, and Cov(E) = Γ =
diag(σ12 , . . . , σp2 ).
17

Theorem 5. Assume Z ∼ N (0, IK ) and E ∼ N (0, Γ) are independent. Let K ≥ 2, p ≥ 2K + 1,
1 ≤ s ≤ 4K/5 and
r
log(K/s)
s
≤ C1 ,
(31)
n
for some universal constant C1 > 0. Then, for all 1 ≤ q ≤ ∞,
)
(
r
log(K/s)
1/q
b A) ≥ C2 s
≥ C3 ,
(32)
inf sup PA Lq (A,
n
b A∈A
A

2
for some positive constants C2 , C3 which only depend on σmin
:= mini∈[p] Γii . The infimum is taken
1/q
b
over all estimators A of A and we use the convention s = 1 for q = +∞.

Since for C = I, the quantity kC −1 k∞,1 becomes 1, we see that the rates of Theorems 4 and
5 coincide up to the logarithmic factors log(K/s) and log(p ∨ n). This shows that our procedure
is not only adaptive, but nearly minimax optimal up to the logarithmic factor provided kC −1 k∞,1
remains bounded.

4.3

Statistical guarantee for Gb and Gbs .

For easy of presentation, and without loss of generality, throughout this section, we continue to
write A for its orthonormal transformation AP that uses the optimal signed permutation matrix
b
P ∈ HK from Theorem 4 to align the columns and signs of A with that of A.
b = K. The latter
We define two criteria to evaluate the estimated clusters Gb on the event K
holds with high probability by Theorem 3. We first define the individual Group False Positive
Proportion (GFPP) and the individual Group False Negative Proportion (GFNP) as,
b a ) :=
GFPP(G

ba |
|(Ga )c ∩ G
,
c
|(Ga ) |

b a ) :=
GFNP(G

b a )c |
|Ga ∩ (G
,
|Ga |

for each a ∈ [K],

(33)

b a )c := [p] \ G
b a , with the convention GFPP(G
b a ) = 0 if |(Ga )c | = ∅.
where (Ga )c := [p] \ Ga and (G
b a . Furthermore,
GFPP and GFNP quantify the misclassification proportion within each group G
with the same convention, we can define the Total False Positive Proportion (TFPP) and Total
b
False Negative Proportion (TFNP) to quantify the overal misclassification proportion of G.
b :=
TFPP(G)

PK

c
b
a=1 |(Ga ) ∩ Ga |
,
PK
c
a=1 |(Ga ) |

b :=
TFNP(G)

Finally, given µ = 5kΩk∞,1 δ 0 with δ 0 specified in (29), we define

PK

b c
a=1 |Ga ∩ (Ga ) |
.
PK
a=1 |Ga |

J2 := {i ∈ J : for any a with Aia 6= 0, |Aia | > (2µ) ∨ (4δ/ν)}.

(34)

(35)

and J3 := J \ (J1 ∪ J2 ). J2 can be viewed as the set where every non-zero entry of Aj· is separated
away from 0 for each j ∈ J2 . The following theorem shows that J2 plays a critical role in quantifying
b and the misclassification proportion of G.
b Let Sb := supp(A).
b
both the support recovery of A
18

Theorem 6. Under the conditions of Theorem 4, with probability greater than 1 − c1 (n ∨ p)−c2 for
some positve constant c1 and c2 defined in (24), we have:
(a)


b ⊆ supp (A) ,
supp AJ2 ⊆ supp(A)


b b) = sign A b .
sign(A
S
S

(b) Let saj = 1{|Aja | =
6 0} and taj = 1{|Aja | ≤ (2µ) ∨ (4δ/ν)}, for each j ∈ J and a ∈ [K].
P
a
a t
b a ) = 0; GFNP(G
b a ) ≤ Pj∈J1 ∪J3 \J1 j .
GFPP(G
a
j∈J sj + |Ia |
(c) Let sj =

Remark 4.

PK

a=1 1{|Aja |

=
6 0} and tj =

b = 0;
TFPP(G)

PK

≤ (2µ) ∨ (4δ/ν)}, for each j ∈ J.
P
j∈J1 ∪J3 tj
b
TFNP(G) ≤ P
.
j∈J sj + |I|

(36)

a=1 1{|Aja |

(37)

1. From our proof of Theorem 6, it is easy to verify that the expression of TFNP in (37) continues to hold for the Direction False Positive Proportion (DFPP) and the Direction False
P
Negative Proportion (DFNP) defined in (40) below with sj replaced by K
a=1 1{Aja < 0} or
P
PK
PK
K
a=1 1{0 < Aja ≤
a=1 1{−(2µ) ∨ (4δ/ν) ≤ Aja < 0} or
a=1 1{Aja > 0}, tj replaced by
+
−
±
(2µ) ∨ (4δ/ν)} and I replaced by I or I , where I := ∪a∈[K] {i ∈ Ia : Aia = ±1}.
b a ) will be small if either ta is small
2. According to display (36), it is easy to see that GFNP(G
j
P
for j ∈ J1 ∪ J3 or |J1 | + |J3 | − |J1a | is dominated by |Ia | + j∈J saj . Moreover, from display
(37), TFNP will be small in the following two cases:
- |J1 | + |J3 | is dominated by |I| + |J2 |;

- tj is small relative to sj , for j ∈ J1 ∪ J3 .
To illustrate this, consider tj ≡ t and sj ≡ s, for each j ∈ J, to simplify the expressions a
bit, and assume |J1 | + |J3 | = α(|I| + |J2 |), for some α ≥ 0. We show in the supplementary
material that



1
(s − 1)|J2 |
b ≤t
TFNP(G)
s+
1+
,
α
|I| + |J2 |

Thus, when either t or α is small, that is, when |J1 | + |J3 | is dominated by |I| + |J2 |, then
TFNP will be small. Note that even when t itself is large but bounded by some constant, TFNP
p
might also be small since s can be close to K which is allowed to grow as O( n/ log p).

3. If J2 = J with µ = 3kC −1 k∞,1 δ, from noting that J2 ⊆ J \ J1 , Remark 2 in Section 4.1 yields
Ib = I. We can choose λ = δ in (15) and µ = 3kC −1 k∞,1 δ in (17), and follow the proof of
Theorems 4 and 6 to arrive at the following conclusions:
b = supp (A) ,
supp(A)

Moreover, we get exact cluster recovery:

19

b = sign(A).
sign(A)

b a ) = GFNP(G
b a ) = 0, for each a ∈ [K].
(a) GFPP(G
b = TFNP(G)
b = 0.
(b) TFPP(G)

b 0 = G0 . Again, all statements hold with probability greater than
This immediately yields G
1 − c1 (n ∨ p)−c2 .

e in which we combine
4. We prove that Theorem 6 also holds for the hard threshold estimator A
j
j
j
j
bb with A
e b. Each row of A
e b is estimated by βea = β̄a 1{|β̄a | > µ} of βa = Aja , a ∈ [K],
b using
A
I
J
J
−1
0
the same µ = 5kC k∞,1 δ as before for the threshold µ. However, we cannot guarantee that
the scaling restriction of condition (i) holds for this estimator.

4.4

Discussion and related work

To the best of our knowledge, optimal estimation of identifiable sparse loading matrices A in model
(1) satisfying (i) - (iii), when both I and K are unknown, and when the entries in X, Z and A are
allowed to have arbitrary signs, has not be considered elsewhere and our results bridge this gap.
There exists, however, a very large body of literature on related problems. We review the most
closely related results below, and explain the differences with our results.
Results regarding the identifiability of A in general latent models, typically without sparsity,
are scattered throughout over more than six decades of literature. They all involve conditions on
both A and C, and there is typically a trade-off between the restrictions on A versus those on C,
as first summarized and proved in Anderson and Rubin [1956], reviewed in Lawley and Maxwell
[1971] and later in Anderson and Amemiya [1988]. We recall them briefly here for the convenience
of the reader.
By far the most commonly used assumption is that the latent factors are uncorrelated, so that
C is either the identity or a diagonal matrix. In this case, it is typically further assumed that the
scaled columns of A are orthogonal, see, for instance, the literature review in Izenman [2008]. An
alternative requirement is that A contain a K × K lower diagonal matrix, see, e.g., Geweke and
Zhou [1996] and, moreover, that the placement of this matrix within A is known, which requires
careful justification [Carvalho et al., 2008], and may be problematic from a practical perspective
[Bhattacharya and Dunson, 2011].
In general, latent factors are correlated, which is our point of view in this work. Then, starting
with Anderson and Rubin [1956], one places on the structure of A constraints that are different than
those made when C is diagonal. The most common of those assumptions involves the existence of
a pure variable set I, similar to our assumption (ii). If I is known, classical results in Anderson
and Rubin [1956] and the proof of our Theorem 2 show that C can be an arbitrary positive definite
matrix. When I is unknown, conditions on the latent factors also need to be imposed. Sufficient
conditions on Z, with provable guarantees for the identification of I, are only known, to the best
of our knowledge, in the NMF literature: the uniqueness of I follows from the uniqueness of the
solution of an appropriate linear program, applied to population quantities, and tailored to matrices
with non-negative entries, see Bittorf et al. [2012]. In contrast, the arguments of Section 2 above
are optimization-free and can be used for matrices that have entries of arbitrary sign. Therefore,

20

we provide a new addition to the literature on pure-variable and loading matrix identification, in
general latent models, and also in the particular case of NMF.
A related, but different, identifiability question regards the covariance matrix Σ of X which,
under (1), can be written as the sum between a rank K matrix and a diagonal matrix:
Σ = ACAT + Γ,

(38)

and Γ := Cov(E) is a diagonal matrix with possibly different entries. In these models, the identifiability question is whether Σ can be decomposed uniquely as the sum between ACAT and Γ.
Answers to this question generated a large amount of literature. We refer the reader to Ledermann
[1937]; Anderson and Rubin [1956]; Shapiro [1982, 1985]; Bekker and ten Berge [1997] for earlier
results, and to Bai and Ng [2002]; Chandrasekaran et al. [2011, 2012]; Candès et al. [2011]; Hsu
et al. [2011], Fan et al. [2013], Wegkamp and Zhao [2016] for more recent works, that also address
the problems of rank estimation and optimal estimation of high dimensional covariance matrices. It
is noteworthy that these works, relative to one another, give different types of sufficient conditions
under which one can separate the low rank matrix ACAT from Γ. However, since we always have
ACAT = (AQ)(QT CQ)(QT AT ), for any orthonormal Q, they do not guarantee the identifiability
of A itself. Conversely, we show in Theorem 2 in Section 2 that under conditions (i) - (iii), C
and Γ are identified, and A is identified up to signed permutations. Therefore, we also identify
uniquely the decomposition of Σ. Our conditions are not always comparable to those employed for
the unique decomposition of Σ, but in special cases they imply them. Although the uniqueness of
the decomposition of Σ is a by-product of our results, we do not pursue the covariance estimation
problem in this work, but we included the above discussion for completeness.
Furthermore, we do not view the problem of estimating the number of factors K as that of
estimating the rank of a matrix. This approach is taken in Bai and Ng [2002], via penalized
least squares, but provided that either C = I or AAT = I and that K is bounded by a fixed
integer. Alternatively, we could adapt the criteria in Bunea et al. [2011]; Bing and Wegkamp
[2018]; Wegkamp and Zhao [2016] to (1) to allow for K → ∞ in the rank estimation problem.
However, proving that such an estimator is consistent would ultimately require an unnecessary
lower bound restriction on the K-th largest eigenvalue of ACAT . In contrast, our Theorem 3 shows
that such conditions can indeed be avoided. We estimate directly the set I and its partition via
LOVE, and as a byproduct K, at a low computational cost of order p2 .
Estimation of A in identifiable factor models is typically based on iterative alternating least
squares procedures or the EM algorithm, see for instance Rubin and Thayer [1982]; Bai and Li
[2012] and the references therein. As discussed in these works, the resulting algorithms are not
suitable for large data sets due to their notoriously slow convergence to a solution that is typically
not the global optimum. Bayesian estimation, see, e.g. Carvalho et al. [2008] and the references
therein, offers an alternative approach which may become computationally very demanding in high
dimensions, requires a likelihood framework, and careful prior specification. Moreover, existing
procedures do not estimate A under our model specifications (i) - (iii), and any adaptation would
still require the challenging estimation of I. Our procedure offers a solution to the computational
problem, as LOVE does not require a likelihood or other prior distributional specifications, is
tailored to our model with unknown I, and has provable low computational complexity.
21

The statistical properties of estimators of A in model (1) (i) - (iii) have not been studied,
and even particular cases of the model have received a very limited amount of attention, from
a theoretical perspective. When I is known and K is fixed, Bai and Li [2012] established the
asymptotic normality of the MLE in a model similar to ours, although the estimator they ultimately
construct is not necessarily the MLE under this model, but rather an appropriate transformation
of the stationary point of a quasi-likelihood for a different factor model. We give the specific details
of their construction in Supplement C.1. If I is unknown, but K is known, and moreover, the
columns of X, A and Z have non-negative entries that sum up to 1, Arora et al. [2013] provide a
practical algorithm for the estimation of A and offer bounds on the `1 matrix norm loss of their
estimator. The extra restrictions on this model are motivated by a specific model, the topic model,
appropriate for vectors with discrete distributions, for instance multinomial. The construction and
analysis of these estimates are not transferable to our general framework, as they depend heavily
on these restrictions. Our results of Section 4.2 bridge this gap in the literature and offer lower and
upper bounds for the performance of estimators of A in model (1) (i) - (iii).
Finally, to the best of our knowledge, overlapping clustering based on model (1) has not been
analyzed. A particular case of this model, corresponding to a matrix A with binary entries, has
been considered in Bunea et al. [2016a,b] for non-overlapping clustering. According to their model,
all p variables are pure variables, as the model assume that Xj = Zk + Ej , for all j ∈ Gk and
k ∈ {1, . . . , K}, {Gk }1≤k≤K form a partition of {1, . . . , p}. When C is positive definite, the nonoverlapping clusters are shown to be identifiable, and the work of Bunea et al. [2016a,b] is devoted
to exact recovery of clusters with minimax optimal cluster separation, a very different problem
than the one considered here.

5

Simulation Studies

In this section, we first discuss our procedure for selecting the tuning parameters, then evaluate the
performance of LOVE based on estimation error and overall clustering misclassification proportion.
In the supplementary materials, we compare LOVE with existing overlapping clustering algorithms
and study the performance of LOVE for the non-overlapping clustering problem.

5.1

Data driven choice of the tuning parameters

Tuning parameter δ. Proposition 3 specifies the theoretical rate of δ, but only up to constants
that depend on the underlying data generating mechanism. We propose below a data-dependent
way to select δ, based on data splitting. Specifically, we split the data set into two independent
b (1) . On the
parts, of equal sizes. On the first set, we calculate the sample covariance matrix Σ
p
second set, we choose a fine grid of values δ` = c` log p/n, with 1 ≤ ` ≤ M , for δ, by varying the
b
proportionality constants c` . For each δ` , we obtain the estimated number of clusters K(`)
and
b
b
b
b
the pure variable set I(`) with its partition I(`). Then we construct the |I(`)| × K(`) submatrix
bb of A,
b and estimate C(`)
b
b
b
A
via formula (11). Finally, we calculate the |I(`)|
× |I(`)|
matrix
I(`)
bb C(`)
b A
bT . In the end, we have constructed a family F = {W1 , . . . , WM } of the fitted
W` = A
I(`)

b
I(`)

22

b
matrices W` , each corresponding to different I(`)
that depend in turn on δ` , for ` ∈ {1, . . . , M }.
Define
1
b (1)
b
,
(39)
CV (I(`))
:= q
b I(`)
b − W`
 ΣI(`)
F-off
b
b
|I(`)|
|I(`)|
−1

where kBkF-off := kB − diag(B)kF denotes the Frobenius norm over the off-diagonal elements
b
of a square matrix B. We choose δ cv as the value δ` that minimizes CV (I(`))
over the grid
` ∈ [M ]. To illustrate how the selection procedure works, we provide an example in Section B of
the supplementary material.
Tuning parameters λ and µ. The tuning parameter λ in the linear program (15) for estimating
Ω = C −1 is specified by λ = 2δ 0 with δ 0 defined in (29). Since δ 0 is proportional to δ, we use λ = c0 δ cv
where c0 is some constant and could be tuned by a cross-validation strategy used in the related
work on the precision matrix estimation, for instance Cai et al. [2011]. More precisely, we randomly
b on the first dataset for each
split the data into two parts. For a given grid of λ, we compute Ω
value in the grid. Then we choose the one which gives the smallest likelihood loss from the second
dataset, where the likelihood loss is defined by
L(Ω, C) = hΩ, Ci − log det(Ω).
From Remark 4 (3) in Section 4.3, when J2 = J, we can choose λ = δ which is the smallest λ
we should consider. Therefore, we set the grid of λ equal to [δ cv , 3δ cv ]. From our simulation, the
selected λ is δ cv in most cases. Hence we recommend to use λ = δ cv and our simulations are based
on this choice.
b estimates C −1 . Our extensive
Recall that µ = c1 kC −1 k∞,1 δ for some constant c1 , and that Ω
b ∞,1 δ cv yields stable performance, with Ω
b solved from
simulations show that the choice of µ = kΩk
(14) and δ cv selected via cross-validation.

5.2

Estimation error and cluster recovery with LOVE

In this section, we study the numerical performance of LOVE in terms of clustering and estimation
accuracy. To the best of our knowledge, there is no comparable algorithm with provable guarantees
developed for our framework, especially if the set I is unknown, as explained in detail in Section
4.4 above, and further re-visited in Section C.1 of the supplementary material.
We generate the data in the following way. We set the number of clusters K to be 20 and
simulate the latent variables Z = (Z1 , . . . , ZK ) from N (0, C). The diagonal elements of C is given
by Cii = 2 + (i − 1)/19 for i = 1, . . . , 20, and the off-diagonal elements are generated as Cij =
(−1)(i+j) 0.3|i−j| (Cii ∧ Cjj ) for any i 6= j. In addition, the error terms E1 , ..., Ep are independently
sampled from N (0, σp2 ), where σp2 itself is sampled from a uniform distribution on [1, 3]. Since the
rows of A corresponding to pure variables in the same cluster are allowed to have different signs, we
consider the following configuration of signs for pure variables in each cluster: (3, 2), (4, 1), (2, 3),
(1, 4) and (5, 0), with the convention that the first number denotes the number of positive pure
variables in that group and the second one denotes the number of negative pure variables. Among
23

the 20 groups, each sign pattern is repeated 4 times. To generate AJ , for any j ∈ J, we randomly
assign the cardinality sj of the support of Aj· to a number in {2, 3, 4, 5}, with equal probability.
Then, we randomly select the support from {1, 2, . . . , K} with cardinality equal to sj . For Ajk
which is nonzero, we set it as Ajk = sign ·(1/sj ) with sign randomly sampled from {−1, 1}. Thus,
we can generate X according to the model X = AZ + E. In the simulation studies, we vary p from
200 to 1000 and n from 300 to 1000. Each simulation is repeated 50 times.
b are not directly comparable, since
Recall that the true allocation matrix A and our estimator A
they may differ by a permutation matrix. To evaluate the performance of our method, we consider
b have the same dimension, we first
the following mapping approach [Wiwie et al., 2015]. If A and A
b kF is minimized.
find the mapping (i.e., the signed permutation matrix P ∈ HK ) such that kA − AP
q = AP
b with A to evaluate the estimation and
Thus, we can compare the permuted estimator A
recovery error. Under this mapping approach, we can evaluate TFPP and TFNP defined in (34).
Moreover, in order to account for the direction sub-groups defined in (21), we can define Direction
False Positive Proportion (DFPP) and Direction False Negative Proportion (DFNP) as follows:
DFPP =

PK

1
b2
a=1 |Ga ∩ Ga |
,
PK
1
a=1 |Ga |

DFNP =

PK

2
b1
a=1 |Ga ∩ Ga |
.
PK
2
a=1 |Ga |

(40)

Figure 1 shows the percentage of exact recovery of number of clusters K, TFPP, TFNP, DFPP
b = K, we can
and DFNP of LOVE. Since the last four measures are well defined only if rank(A)
compute them when the number of clusters is correctly identified. We can see that the proposed
method correctly selects K and as long as the number of clusters is correctly selected, TFPP, TFNP,
DFPP and DFNP of our method are very close to 0, which implies that the sign and sparsity pattern
b as measured by the matrix
of A can be correctly recovered. We present the estimation error of A
√
`1 norm scaled by pK and the Frobenius norm scaled by pK in Table 1.
As expected, the estimation error decreases when the sample size increases from 300 to 1000,
which is in line with our theoretical results. The simulations are conducted on an macOS Sierra
system version 10.12.6 with 2.2 GHz Intel Core i7 CPU and 16 GB memory. Even with p = 1000
and n = 1000, the computing time of our method for each simulation is around 1 minute.
Moreover, we evaluated the performance of the LOVE procedure for K varying in a wide range,
from 3 to 30, and when AJ contains many very small entries. The results are consistent with what
we observed in this section and deliver the same message. The GFPP and GFNP are similar as
e defined in Remark 4
TFPP and TFNP and the performance of the hard thresholding estimator A,
b To save space, we have omitted those results.
of Section 4.3, is similar to A.
We also compared the performance of LOVE with other off-the-shelf algorithms for overlapping
clustering, and tested LOVE for non-overlapping clustering. We included these results in sections
C.2 and C.3 of the supplementary material.

6

Application

To benchmark LOVE, we used a publicly available RNA-seq dataset of 285 blood platelet samples
from patients with different malignant tumors [Best et al., 2015]. We extracted a small subset of 500
24

b as measured by the matrix `1 norm (`1 ) (divided by
Table 1: The average estimation error of A
√
pK) and the Frobenius norm (`2 ) (divided by pK). Numbers in parentheses are the simulation
standard errors.

p
200
400
600
800
1000

n = 300

n = 500

n = 700

`1

`2

`1

`2

`1

`2

`1

`2

0.018
(0.001)
0.026
(0.002)
0.029
(0.002)
0.031
(0.002)
0.032
(0.002)

0.062
(0.005)
0.075
(0.007)
0.079
(0.006)
0.083
(0.006)
0.083
(0.006)

0.015
(0.001)
0.023
(0.001)
0.025
(0.001)
0.026
(0.001)
0.027
(0.001)

0.053
(0.003)
0.064
(0.003)
0.067
(0.003)
0.068
(0.004)
0.069
(0.003)

0.013
(0.001)
0.021
(0.001)
0.023
(0.001)
0.024
(0.001)
0.025
(0.001)

0.048
(0.008)
0.059
(0.006)
0.063
(0.003)
0.064
(0.004)
0.065
(0.004)

0.012
(0.001)
0.018
(0.001)
0.020
(0.001)
0.022
(0.001)
0.022
(0.001)

0.041
(0.002)
0.051
(0.003)
0.055
(0.003)
0.057
(0.004)
0.057
(0.004)

0.8
0.6
0.4
0.2
0.0

0.0

0.2

cluster
TFPP
TFNP
DFPP
DFNP

recovery and error rates

0.4

0.6

0.8

1.0

p = 1000

1.0

p = 400

recovery and error rates

n = 1000

300 400 500 600 700 800 900

300 400 500 600 700 800 900

sample size

sample size

Figure 1: Percentage of exact recovery of number of clusters K (cluster), total false positive proportion (TFPP), total false negative proportion (TFNP), direction false positive proportion (DFPP)
and direction false negative proportion (DFNP) for LOVE.
Ensembl genes to test the method. The goal of the benchmarking was to test whether (i) clusters
corresponded to biological knowledge, specifically Gene Ontology (GO) functional annotation of the
genes [Ashburner et al., 2000], (ii) overlapping clusters corresponded to pleiotropic gene function.
LOVE produced twelve overlapping clusters (Table 2) which aligned well with a-priori expectation.

25

Table 2 lists the number of pure genes and the total number of genes in twelve overlapping clusters.
Figure 2 shows that each cluster overlaps with the other and also gives us a clear picture on how two
clusters possibly overlap. For example, 18 genes belong to both cluster 3 and cluster 11, whereas
cluster 2 and cluster 3 have only one common gene. The genes with the same GO biological
process, molecular function or cellular component terms tended to be assigned to the same cluster.
For example, ENSG00000273906 and ENSG00000273328 are both RNA genes. They were both
assigned to the same cluster (cluster 6, Figure 2). However, they were also assigned to other
clusters, suggesting they have pleiotropic functions. This suggests that the latent variables used
for clustering are likely to have biological significance and can potentially be used for functional
discovery for genes with under-explored functions. We found 308 genes with zero expression across
all samples. None of them were assigned to any of the 12 estimated clusters, as desired. Indeed,
our model not only allows for the existence of pure noise variables Xj = Ej , but variables with
structural zero values as well, as Γjj = Var(Ej ) = 0 is permitted. Formally we place them in the
pure noise cluster G0 , for further scientific scrutiny.
Table 2: Number of pure genes and total number of genes in each group.
G1

G2

G3

G4

G5

G6

G7

G8

G9

G10

G11

G12

2
58

2
35

2
67

4
105

2
80

10
104

2
28

2
43

2
44

4
74

2
94

15
108

Number of pure genes
Total number of genes
14

8

G67
G5

G3

10

11

6

9

7

6

15

Group size

13

28

18

35

10
5

43

1
6

G4

9

8

44

7
7

3

8
4

8

4

G2

11 G8
2
14 7
4

4

11

3

102 4

8

7

14
9

G7

4

1

80
94

3 10

3
1

74

6

14

7

2

67

9
5

4

58

G11

16

6

7

3

104
105

1
3

3

3

9

G9
G10 6

108

7
4

G1

5

G12

Figure 2: Left panel: Number of genes overlapped in different groups. The nodes represent 12
groups with the same labels and sizes as those in Table 2. The number shown on the edge between
two nodes represents the number of genes shared by the two groups, which corresponds to the width
of that edge. Right panel: Illustration of three genes ENSG00000273906, ENSG00000273328 and
ENSG00000273113 and their allocation matrix relative to 12 groups. For instance, the jth gene
bj6 = 0.04, A
bj9 = 0.37, A
bj11 = −0.02.
ENSG00000273906 belongs to groups 6, 9 and 11 with A
26

Acknowledgements
We are grateful to Jishnu Das for help with the interpretation of our data analysis results. Bunea
and Wegkamp are supported in part by NSF grant DMS 1712709. Bing was supported in part by
NSF grant DMS-1407600.

Supplementary Material
The supplementary materials contain the proofs and additional numerical results.

References
Anderson, T. (2003). An Introduction to Multivariate Statistical Analysis. Wiley Series in Probability and Statistics, Wiley.
Anderson, T. W. and Amemiya, Y. (1988). The asymptotic normal distribution of estimators
in factor analysis under general conditions. Ann. Statist. 16 759–771.
Anderson, T. W. and Rubin, H. (1956). Statistical inference in factor analysis. In Proceedings
of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 5: Contributions to Econometrics, Industrial Research, and Psychometry. University of California Press,
Berkeley, Calif.
Arora, S., Ge, R., Halpern, Y., Mimno, D. M., Moitra, A., Sontag, D., Wu, Y. and Zhu,
M. (2013). A practical algorithm for topic modeling with provable guarantees. In ICML (2).
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M.,
Davis, A. P., Dolinski, K., Dwight, S. S. and Eppig, J. T. (2000). Gene ontology: tool for
the unification of biology. Nature genetics 25 25–29.
Bai, J. and Li, K. (2012). Statistical analysis of factor models of high dimension. Ann. Statist.
40 436–465.
Bai, J. and Ng, S. (2002). Determining the number of factors in approximate factor models.
Econometrica 70 191–221.
Bekker, P. A. and ten Berge, J. M. (1997). Generic global identification in factor analysis.
Linear Algebra and its Applications 264 255 – 263. Sixth Special Issue on Linear Algebra and
Statistics.
Belloni, A., Rosenbaum, M. and Tsybakov, A. B. (2017). Linear and conic programming
estimators in high dimensional errors-in-variables models. Journal of the Royal Statistical Society:
Series B (Statistical Methodology) 79 939–956.
Best, M. G., Sol, N., Kooi, I., Tannous, J., Westerman, B. A., Rustenburg, F.,
Schellen, P., Verschueren, H., Post, E., Koster, J. et al. (2015). Rna-seq of tumoreducated platelets enables blood-based pan-cancer, multiclass, and molecular pathway cancer
diagnostics. Cancer cell 28 666–676.
Bezdek, J. C. (2013). Pattern recognition with fuzzy objective function algorithms. Springer
Science & Business Media.

27

Bhattacharya, A. and Dunson, D. (2011). Sparse bayesian infinite factor models. Biometrika
98 291–306.
Bien, J., Bunea, F. and Xiao, L. (2016). Convex banding of the covariance matrix. Journal of
the American Statistical Association 111 834–845.
Bing, X. and Wegkamp, M. H. (2018). Adaptive estimation of the rank of the coefficient matrix
in high dimensional multivariate response regression models. ArXiv .
Bittorf, V., Recht, B., Re, C. and Tropp, J. A. (2012). Factoring nonnegative matrices with
linear programs. arXiv:1206.1270 .
Bollen, K. A. (1989). Structural Equations with Latent Variables. Wiley.
Bunea, F., Giraud, C. and Luo, X. (2016a). Minimax optimal variable clustering in g-models
via cord. arXiv preprint arXiv:1508.01939 .
Bunea, F., Giraud, C., Royer, M. and Verzelen, N. (2016b). Pecok: a convex optimization
approach to variable clustering. arXiv preprint arXiv:1606.05100 .
Bunea, F., She, Y. and Wegkamp, M. H. (2011). Optimal selection of reduced rank estimators
of high-dimensional matrices. Ann. Statist. 39 1282–1309.
Cai, T., Liu, W. and Luo, X. (2011). A constrained 1 minimization approach to sparse precision
matrix estimation. Journal of the American Statistical Association 106 594–607.
Cai, T. T., Liu, W. and Zhou, H. H. (2016). Estimating sparse precision matrix: Optimal rates
of convergence and adaptive estimation. Annals of Statistics 44 455–488.
Candès, E. J., Li, X., Ma, Y. and Wright, J. (2011). Robust principal component analysis?
Journal of the ACM (JACM) 58 11.
Carvalho, C. M., Chang, J., Lucas, J. E., Nevins, J. R., Wang, Q. and West, M. (2008).
High-dimensional sparse factor modeling: Applications in gene expression genomics. Journal of
the American Statistical Association 103 14381456.
Chandrasekaran, V., Parrilo, P. A. and Willsky, A. S. (2012). Latent variable graphical
model selection via convex optimization. Ann. Statist. 40 1935–1967.
Chandrasekaran, V., Sanghavi, S., Parrilo, P. A. and Willsky, A. S. (2011). Ranksparsity incoherence for matrix decomposition. SIAM J. Optim. 21 572–596.
Craddock, R. C., James, G. A., Holtzheimer, P. E., Hu, X. P. and Mayberg, H. S.
(2012). A whole brain fmri atlas generated via spatially constrained spectral clustering. Human
brain mapping 33 1914–1928.
Craddock, R. C., Jbabdi, S., Yan, C.-G., Vogelstein, J. T., Castellanos, F. X., Di Martino, A., Kelly, C., Heberlein, K., Colcombe, S. and Milham, M. P. (2013). Imaging
human connectomes at the macroscale. Nature methods 10 524–539.
Donoho, D. and Stodden, V. (2004). When does non-negative matrix factorization give a correct
decomposition into parts? In Advances in Neural Information Processing Systems 16 (S. Thrun,
L. K. Saul and P. B. Schölkopf, eds.). MIT Press, 1141–1148.
Everitt, B. (1984). An Introduction to Latent Variable Models. Monographs on Statistics and
Applied Probability. Springer.
Fan, J., Liao, Y. and Mincheva, M. (2013). Large covariance estimation by thresholding prin-

28

cipal orthogonal complements. Journal of the Royal Statistical Society: Series B (Statistical
Methodology) 75 603–680.
Friedman, J., Hastie, T. and Tibshirani, R. (2008). Sparse inverse covariance estimation with
the graphical lasso. Biostatistics 9 432–441.
Gautier, E. and Tsybakov, A. B. (2011). High-dimensional instrumental variables regression
and confidence sets. arXiv preprint arXiv:1105.2454v4 1105.2454.
Geweke, J. and Zhou, G. (1996). Measuring the pricing error of the arbitrage pricing theory.
The review of financial studies 9 557–587.
Hsu, D., Kakade, S. M. and Zhang, T. (2011). Robust matrix decomposition with sparse
corruptions. IEEE Trans. Inform. Theory 57 7221–7234.
Izenman, A. J. (2008). Modern Multivariate Statistical Techniques: Regression, Classification,
and Manifold Learning. Series: Springer Texts in Statistics.
Jiang, D., Tang, C. and Zhang, A. (2004). Cluster analysis for gene expression data: A survey.
IEEE Transactions on knowledge and data engineering 16 1370–1386.
Koopmans, T. C. and Reiersol, O. (1950). The identification of structural characteristics. Ann.
Math. Statist. 21 165–181.
Krishnapuram, R., Joshi, A., Nasraoui, O. and Yi, L. (2001). Low-complexity fuzzy relational
clustering algorithms for web mining. IEEE transactions on Fuzzy Systems 9 595–607.
Lawley, D. N. and Maxwell, A. E. (1971). Factor analysis as a statistical method. 2nd ed.
American Elsevier Publishing Co., Inc., New York.
Ledermann, W. (1937). On the rank of the reduced correlational matrix in multiple-factor analysis. Psychometrika 2 85–93.
McDonald, R. P. (1999). Test theory: a unified treatment. Taylor and Francis.
Meinshausen, N. and Bühlmann, P. (2006). High-dimensional graphs and variable selection
with the lasso. Annals of Statistics 1436–1462.
Rigollet, P. and Tsybakov, A. (2011). Exponential screening and optimal rates of sparse
estimation. Ann. Statist. 39 731–771.
Rubin, D. B. and Thayer, D. T. (1982). Em algorithms for ml factor analysis. Psychometrika
47 69–76.
Shapiro, A. (1982). Rank-reducibility of a symmetric matrix and sampling theory of minimum
trace factor analysis. Psychometrika 47 187–199.
Shapiro, A. (1985). Identifiability of factor analysis: some results and open problems. Linear
Algebra and Its Applications 70 1–7.
Tsybakov, A. B. (2009). Introduction to Nonparametric Estimation. Springer, New York.
Wegkamp, M. and Zhao, Y. (2016). Adaptive estimation of the copula correlation matrix for
semiparametric elliptical copulas. Bernoulli 22 1184–1226.
Wiwie, C., Baumbach, J. and Röttger, R. (2015). Comparing the performance of biomedical
clustering methods. Nature methods 12 1033–1038.
Yuan, M. and Lin, Y. (2007). Model selection and estimation in the gaussian graphical model.
Biometrika 94 19–35.

29

A
A.1

Appendix
Proofs of the results from Section 2

We begin by stating and proving two lemmata that are crucial for the main results of this section.
All results are proved under the condition that model 1 and (i) - (iii) hold.
Lemma 7. For any a ∈ [K] and i ∈ Ia , we have
(a) |Σij | = Caa for all j ∈ Ia ,
(b) |Σij | < Caa for all j 6∈ Ia .
Proof. For given i ∈ [p], we define the set s(i) := {1 ≤ a ≤ K : Aia 6= 0}. For any i ∈ Ia and j 6= i,
we have
|Σij | =
=

X

Aia

a∈s(i)

X

X

Ajb Cab

b∈s(j)



Ajb Cab from the definition of Ia

b∈s(j)

≤

X

b∈s(j)

|Ajb | · max |Cab |
b∈s(i)

≤ Caa using conditions (i) and (iii).
Furthermore, using conditions (i) and (iii), we observe that we have equality in the above display
for j ∈ Ia , and strict inequality for j 6∈ Ia , which proves the lemma.
Lemma 8. We have
(a) Si ∩ I 6= ∅, for any i ∈ [p],
(b) Si ∪ {i} = Ia and Mi = Caa , for any i ∈ Ia and a ∈ [K],
where Mi and Si are defined in (4) and (5), respectively.
Proof. Lemma 7 implies that, for any i ∈ Ia , Mi = Caa and Si = Ia \ {i}, which proves part (b).
From the result of part (b), it remains to show Si ∩ I 6= ∅ for any i ∈
/ I. Let i ∈
/ I be fixed. We
have
!
X
X
(41)
Ajb
Aia Cab
Mi = max |Σij | = max
j6=i

≤ max max

j6=i b∈s(j)

j6=i

X

b∈s(j)

a∈s(i)

Aia Cab = max
j6=i

a∈s(i)

P

X

Aia Cab∗

a∈s(i)

for some b∗ ∈ [K]. A direct computation yields |Σij | = | a∈s(i) Aia Cab∗ | for any j ∈ Ib∗ , that is,
the maximum Mi of |Σij | is achieved at all j ∈ Ib∗ . Since Ib∗ 6= ∅ by condition (ii), this completes
the proof of claim (a).
30

Proof of Theorem 1. We have all the necessary ingredients to proceed with the proof of the
main result of this section.
Proof of (a). We first show the sufficiency part. Consider any i ∈ [p] with Mi = Mj for all j ∈ Si .
Part (a) of Lemma 8 states that there exists a j ∈ Ia ∩ Si for some a ∈ [K]. For this j ∈ Ia , we
have Mj = Caa from part (b) of Lemma 8. Invoking our premise Mj = Mi as j ∈ Si , we conclude
that Mi = Caa , that is, maxk6=i |Σik | = Caa . By Lemma 7, the maximum is achieved for any pair
i, k ∈ Ia . However, if i 6∈ Ia , we have that |Σik | < Caa for all k 6= i. Hence i ∈ Ia and this concludes
the proof of the sufficiency part.
It remains to prove the necessity part. Let i ∈ Ia for some a ∈ [K] and j ∈ Si . Lemma
8 implies that j ∈ Ia and Mi = Caa . Since j ∈ Si , we have |Σij | = Mi = Caa , while j ∈ Ia
yields |Σjk | ≤ Caa for all k 6= j, and |Σjk | = Caa for k ∈ Ia , as a result of Lemma 7. Hence,
Mj = maxk6=j |Σjk | = Caa = Mi for any j ∈ Si , which proves our claim.
Proof of (b). We start with the following constructive approach. Let N = [p] be the set of all
variable indices and O = ∅. Let Mi and Si be defined in (4) and (5), respectively.
(1) Choose i ∈ N and calculate Si and Mi .
(a) If Mi = Mj , for all j ∈ Si , set I (i) := Si ∪ {i}, O = O ∪ {i} and N = N \ I (i) .

(b) Otherwise, replace N by N \ {i}.
(2) Repeat step (1) until N = ∅.

We show that {I (i) : i ∈ O} = I. Let i ∈ O be arbitrary fixed. By (a), we have i ∈ I. Thus, there
exists a ∈ [K] such that i ∈ Ia . By Lemma 8, i ∈ Ia implies Ia = Si ∪ {i} = I (i) . On the other
hand, let a ∈ [K] be arbitrary fixed. By condition (ii), there exists at least one j ∈ Ia . Once again,
by part (b) of Lemma 8, if j ∈ Ia , then Sj ∪ {j} = Ia , that is, I (j) = Ia .

Proof of Theorem 2
Theorem 1 shows that Σ uniquely defines I and its partition I, up to permutation of labels. Given
I and its partition I = {I1 , . . . , IK }, for any i ∈ I, there exists a unique 1 ≤ a ≤ K such that
i ∈ Ia . Then we set |Ai· | = ea , the canonical basis vector in RK that contains 1 in position a and is
zero otherwise. Thus, the |I| × K matrix AI with rows Ai· is uniquely defined up to multiplication
with a signed permutation matrix P .
We show below that AJ is also identifiable up to a signed permutation matrix. We begin by
observing that, for each i ∈ Ik , for some k ∈ [K], and any j ∈ J, Model 1 implies
X X
X
Σij =
Aia Ajb Ckb = Aik
Ajb Ckb
a∈s(i) b∈s(j)

and since A2ik = 1, we obtain

b∈s(j)

T
Aik Σij = Ck·
Aj·

31

and, after averaging over all i ∈ Ik ,
T
Ck·
Aj· =

1 X
Aik Σij .
|Ik |
i∈Ik

Repeating this for every k ∈ [K], we obtain the formula

T
X
X
1
1
CAj· = 
Ai1 Σij , . . . ,
AiK Σij  := θj .
|I1 |
|IK |
i∈I1

i∈IK

The covariance matrix C can be uniquely constructed from Σ via
1
|Ia |(|Ia | − 1)

Caa =
for any a ∈ [K], and

1
|Ia ||Ib |

Cab =

X

X

i,j∈Ia ,i6=j

|Σij |

Aia Ajb Σij

i∈Ia ,j∈Ib

for a, b ∈ [K] with a 6= b. Notice that mina∈[K] |Ia | ≥ 2, which is part of our model requirement
(ii), is needed for the construction of Caa . Since the covariance matrix C is assumed to be positive
definite, Aj· = C −1 θj , for each j ∈ J, which shows that AJ can be determined uniquely from Σ up
to a signed permutation. Therefore, AJ is identifiable which concludes the proof.

A.2

Proofs of the results from Section 4.1

The proof of Theorem 3 will repeatedly use Lemma 9, stated and proved below. Let
ci := max |Σ
b ij |.
M
j∈[p]\{i}

(42)

Lemma 9. Under the conditions in Theorem 3, for any i ∈ Ia with some a ∈ [K], the following
inequalities hold on the event E:
b ij | − |Σ
b ik | ≤ 2δ,
|Σ
b ij | − |Σ
b ik | > 2δ,
|Σ
b ij | − |Σ
b ik | < 2δ,
|Σ

For any i ∈ J1a , we have

for all j, k ∈ Ia \ {i} and j 6= k;

(43)

for all j ∈ Ia \ {i}, k ∈
/ (Ia ∪ J1a );

(44)

for all j ∈ J1a and k ∈ Ia \ {i}.

ci − |Σ
b ij | ≤ 2δ,
M

for any j ∈ Ia .

(45)

(46)

Proof of Lemma 9. For the entire proof, we work on the event E defined in (23). To prove (43), we
observe that, for any i, j, k ∈ Ia , Σij = Σik = Caa by Lemma 7, whence
b ij | − |Σ
b ik |
|Σ

≤

|Σij | − |Σik | + 2δ = 2δ.
32

To prove (44), we first observe that, for any j ∈ Ia , |Σij | = Caa by Lemma 7, whence
E

b ij | ≥ Caa − δ.
|Σ

Next, we notice that, for any ` ∈ [p],
|Σi` | =

K
X

A`b Cab = A`a Caa +

b=1

(iii)

X

(47)

A`b Cab

b6=a

≤ |A`a |Caa + (1 − |A`a |)(Caa − ν) = Caa − (1 − |A`a |)ν.

(48)

For any j ∈ Ia and k ∈ [p] \ (Ia ∪ J1a ), the definition of J1 implies |Aka | ≤ 4δ/ν, hence
(48)

E

b ik | ≤ |Σik | + δ ≤ Caa − (1 − |Aka |)ν + δ ≤ Caa − ν + 5δ,
|Σ

so that

E

b ij | − |Σ
b ik | ≥ |Σij | − δ − |Σ
b ik | ≥ |Σij | − Caa + ν − 6δ > 2δ,
|Σ

by using ν > 8δ · (kCk∞ /ν) ≥ 8δ. To prove (45), observe that, for any j ∈ J1a and k ∈ Ia \ {i},
(48)

E

b ij | ≤ Caa − (1 − |Aja |)ν + δ < Caa + δ = |Σik | + δ ≤ |Σ
b ik | + 2δ.
|Σ

So far, we have proved (43) - (45) and it remains to show (46). For any i ∈ J1a , we have, for some
c ∈ [K],
E

ci ≤ max |Σik | + δ (41)
M
=
k∈[p]\i

K
X
b=1

(∗)

Aib Cbc + δ ≤

K
X
b=1

E

b ij | + 2δ.
Aib Cba + δ = |Σij | + δ ≤ |Σ

It remains to show that inequality (∗) holds, for any c 6= a. On the one hand, we have
K
X
b=1

(iii)

Aib Cbc ≤ |Aia ||Cac | + (1 − |Aia |)Ccc ≤ |Aia |(Caa − ν) + (1 − |Aia |)Ccc ,

while on the other hand, we find
K
X
b=1

(iii)

Aib Cab ≥ |Aia ||Caa | − (1 − |Aia |)(Caa − ν).

Combining the preceding two display yields
K
X
b=1

Aib Cab −

K
X
b=1

Aib Cbc ≥ ν − (1 − |Aia |)(Caa + Ccc ).

The term on the right is positive, since condition (28) guarantees that
4δ
(Caa + Ccc ) ≥ (1 − |Aia |)(Caa + Ccc ),
ν
where the last inequality is due to the definition of J1 . This concludes the proof.
ν>

33

Lemma 9 remains valid under the conditions of Remark 2 in which case J1 = ∅ and we only
need ν > 4δ to prove (44).
Proof of Theorem 3. We work on the event E throughout the proof. Without loss of generality,
we assume that the label permutation π is the identity. We start by pointing out that the following
three claims are sufficient to prove (a) - (c). Let Ib(i) be defined in step 4 of Algorithm 1.
(1) For any i ∈ J \ J1 , we have P ure(i) = F alse.

(2) For any i ∈ Ia and a ∈ [K], we have P ure(i) = T rue, Ia ⊆ Ib(i) and Ib(i) \ Ia ⊆ J1a .

(3) For any i ∈ J1a and a ∈ [K], we have Ia ⊆ Ib(i) .

If we can prove these claims, then (1) implies that none of variables in J \ J1 will be selected in
any set of Ib via i ∈ J \ J1 . (2) implies that for any a ∈ [K], there exists Iba such that Ia ⊆ Iba and
Iba \ Ia ⊆ J1a . Moreover, this together with Merge in Algorithm 1 prevents Iba from selecting any
variable from [p] \ (Ia ∪ J1a ). Finally, (3) guarantees that none of pure variables will be excluded
b = K and Ib = {Ib1 , . . . , IbK } is the desired partition.
by any i ∈ J1 in the Merge step. Thus, K
Therefore, in the following we proceed to prove (1) - (3).
To prove (1), let i ∈ J \ J1 be fixed. We first prove that P ure(i) = F alse when Ib(i) ∩ I 6= ∅. It
suffices to show that, there exists j ∈ Ib(i) such that the following does not hold
cj − |Σ
b ij | ≤ 2δ.
M

(49)

Let Ib(i) ∩ I 6= ∅, so there exists j ∈ Ib ∩ Ib(i) for some b ∈ [K]. For such j, we have |Σij | =
P
| K
a=1 Aia Cab | and
E

b ij | ≤
|Σ

K
X
a=1

(iii)

Aia Cab + δ ≤ |Aib |Cbb + (1 − |Aib |)(Cbb − ν) + δ < Cbb − 3δ,

(50)

using the definition of J1 to justify the last inequality. On the other hand, since j ∈ Ib , part (b) of
Lemma 8 implies
E
cj = max |Σ
b jk | ≥ max |Σjk | − δ = Cbb − δ.
M
(51)
k∈[p]\{i}

k∈[p]\{i}

cj − |Σ
b ij | > 2δ. This shows that for any i ∈ J \ J1 , if Ib(i) ∩ I 6= ∅,
Combining (50) with (51) gives M
then P ure(i) = F alse. Therefore, to complete the proof of (1), we show Ib(i) ∩ I = ∅ is impossible
when i ∈ J \ J1 under our assumptions. If Ib(i) ∩ I = ∅, then there exists some j ∈ J ∩ Ib(i) and
|Σij | =

K X
K
X
b=1 a=1

Aia Ajb Cab ≤ max

1≤b≤K

K
X

Aia Cab =

a=1

K
X
a=1

Aia Cab∗ = |Σik |

for some b∗ ∈ [K] and any k ∈ Ib∗ (the set Ib∗ is non-empty by condition (ii)). Therefore,
E

b ij | − |Σ
b ik | ≤ |Σij | − |Σik | + 2δ
|Σ
34

≤

2δ

However, since Ib(i) ∩ I = ∅ and k ∈ Ib∗ , we know k ∈
/ Ib(i) , which implies
b ij | − |Σ
b ik | > 2δ,
|Σ

from Step 4 of Algorithm 1. The last two displays contradict each other, and we conclude that, for
any i ∈ J \ J1 , Ib(i) ∩ I 6= ∅.
To prove (2), let i ∈ Ia be arbitrarily fixed with some a ∈ [K]. We first show that P ure(i) =
T rue. From steps 7 - 8 of Algorithm 1, it suffices to show that, for any j ∈ Ib(i) , (49) holds. From
(44) in Lemma 9, given Step 4 of Algorithm 1, we know that, for any j ∈ Ib(i) , j ∈ Ia ∪ J1a . Thus, we
cj is achieved by
write Ib(i) = (Ib(i) ∩ Ia ) ∪ (Ib(i) ∩ J1a ). For any j ∈ Ib(i) ∩ Ia , by the same reasoning, M
a
some element in either Ia or J1 . For both cases, since i, j ∈ Ia and i 6= j, (43) and (45) in Lemma
9 guarantee that (49) holds. On the other hand, for any j ∈ Ib(i) ∩ J1a , (46) in Lemma 9 implies
that (49) still holds. Thus, we have shown that, for any i ∈ Ia , P ure(i) = T rue. To show Ia ⊆ Ib(i) ,
ci can only be achieved by indices in Ia ∪ J a . In both cases,
let any j ∈ Ia \ {i} and observe that M
1
(i)
(43) and (45) imply j ∈ Ib . Thus, Ia ⊆ Ib(i) . Finally, Ib(i) \ Ia ⊆ J1a follows immediately from (44).
We conclude the proof by noting that (3) immediately follows from (46).

A.3

Proofs of the results from Section 4.2

We divide the proof of Theorem 4 into three steps:
bI aligns
Step 1. We show that there exists a signed permutation Pb such that the columns of A
with those of AI in terms of label and sign, as detailed in Lemma 10;
bb − Ā b;
Step 2. We write Ā = APb, and prove first the error bounds for A
I
I
b
b
b
Step 3. We prove the error bounds for A − Ā = A − AP , with the same Pb, which further implies
b and A.
that Pb aligns the columns of A
Lemma 10. Under conditions of Theorem 4, there exists a signed permutation matrix Q such that
bia ) for any i ∈ Iba with each a ∈ [K].
Ā = AQ satisfies that sign(Āia ) = sign(A

b = K, I ⊆ Ib ⊆ I ∪ J1 and Iπ(a) ⊆ Iba ⊆ Iπ(a) ∪ J π(a) ,
Proof of Lemma 10. Theorem 3 guarantees K
1
with high probability, for any a ∈ [K] and some label permutation π. Let us write Q = Q1 Q2 , with
bI ,
the unsigned permutation matrix Q1 which relabels the columns of AI according to those of A
and with Q2 = diag(q1 , . . . , qK ) with qa ∈ {+1, −1} for each a ∈ [K].
q = AQ1 , we proceed to show that, for each a ∈ [K], sign(A
qia ) · qa holds
bia ) = sign(A
Denoting A
π(a)
b
b
for any i ∈ Ia , in which case each qa can be uniquely constructed. Since Ia ⊆ Iπ(a) ∪ J1 , it suffices
to prove that, for any a ∈ [K],
bja )
bia )
sign(A
sign(A
=
,
qia )
qja )
sign(A
sign(A

π(a)

for any i, j ∈ Iπ(a) or i, j ∈ J1

with i 6= j.

(52)

bI , for any i, j ∈ Iπ(a) or i, j ∈ J π(a) , we
From the definition of AI and the way we construct A
1
consider the following two cases:
qia ) = sign(A
qja ). To show A
bia = A
bja , from (6),
If sign(Aiπ(a) ) = sign(Ajπ(a) ), this implies sign(A
1
2
b
b
b
we need to show i, j ∈ Ia or i, j ∈ Ia which is equivalent to show Σij > 0. For any i, j ∈ Iπ(a) or
35

P
π(a)
i, j ∈ J1
with i 6= j, display (26) gives |Akπ(a) | ≥ 1 − 4δ/ν and b6=π(a) |Akb | ≤ 4δ/ν, for k = i, j.
Thus, using sign(Aiπ(a) ) = sign(Ajπ(a) ), we have
X
X
X
Σij = Aiπ(a) Ajπ(a) Cπ(a)π(a) + Aiπ(a)
Ajc Cπ(a)c + Ajπ(a)
Aib Ajc Cbc
Aib Cπ(a)b +
c6=a

b6=a

b,c6=π(a)

(iii)

≥ Aiπ(a) Ajπ(a) Cπ(a)π(a) − |Aiπ(a) |(1 − |Ajπ(a) |)(Cπ(a)π(a) − ν)
X
Aib Ajc Cbc
− |Ajπ(a) |(1 − |Aiπ(a) |)(Cπ(a)π(a) − ν) −


b,c6=π(a)

2



8δ
16δ 2
4δ
4δ
≥ 1−
Cπ(a)π(a) −
Cπ(a)π(a) − 2 Cb∗ b∗ + 8δ,
· 1−
ν
ν
ν
ν


2
2
16δ
16δ
48δ
≥ 1+ 2 −
Cπ(a)π(a) − 2 Cb∗ b∗ + 8δ.
ν
ν
ν

for some b∗ 6= π(a)

b ij ≥ Σij − δ > 3δ > 0.
Since (28) implies 8δCb∗ b∗ < ν 2 and ν > 8δ, on the event E, we have Σ
qia ) 6= sign(A
qja ). Similarly, to show A
bia 6= A
bja ,
If sign(Aiπ(a) ) 6= sign(Ajπ(a) ), this gives sign(A
b ij < 0. Using the same arguments yields
we prove Σ
E

b kl ≤ Σkl + δ < −3δ < 0.
Σ

Therefore, given Ib = {Iba }a∈[K] , we can construct the signed permutation Pb = Q which alligns the
bI .
columns of AI with those of A

For ease of notation and without loss of generality, we make the blanket assumption that the
signed permutation Pb is the identity so that Ā = A for the remainder of the proof. We note that
the signed permutation Pb will be the same when estimating each row Aj· for j ∈ J.

bb and parts (a) - (c) in Theorem 3, we can write, for
Proof of step 2: From the construction of A
I
each a ∈ [K], Iba = Ia ∪ La with La := Iba ∩ J1a . For any i ∈ Iba , the definitions of I and J1a imply
bia ), we have
|Aia | ≥ 1 − 4δ/ν. Since Lemma 10 guarantees that sign(Aia ) = sign(A
bb − A bk∞ = max kA
bi. − Ai. k∞ ≤ 4 δ.
kA
I
I
ν
b
i∈I

b we have
Let si = kAi· k0 for i ∈ [p]. Then, for any i ∈ I,
bj· − Āj·
A

q

≤

4 1/q
s δ,
ν i

1 ≤ q ≤ ∞.

For Step 3 of the proof of Theorem 4, we will make use of the results of Lemmas 11 and 12, stated
here first and proved at the end of this section, in order to preserve the flow of the presentation.
Lemma 11. Under the conditions of Theorem 4, on the event E, we have

where δ 0 is given in (29).

b − Ck∞ ≤ 2δ 0 ,
kC

max kθbj − θj k∞ ≤ δ 0 ,
j∈Jb

36

(53)

Lemma 12. Under the conditions of Theorem 4, on the event E, we have βaj = 0 implies βbaj = 0,
b
for any j ∈ Jb and a ∈ [K].
b recall that β j = C −1 θj = Ωθj since C is invertible. Also recall
Proof of Step 3. For each j ∈ J,
b θbj . We first show kβ̄ j − β j k∞ ≤ 5kΩk∞,1 δ 0 . For notational convenience, we remove all
that β̄ j = Ω
the super indices. From Lemma 11, the following event
(
)
b − Ck∞ ≤ 2δ 0 , max kθbj − θj k∞ ≤ δ 0 ,
E 0 = kC
j∈Jb

is implied by the event E = Eδ . On the event E 0 , the true Ω := C −1 satisfies the constraint since
b − Ik∞ = kΩ(C
b − C)k∞ ≤ kC
b − Ck∞ kΩk∞,1 ≤ 2δ 0 kΩk∞,1 .
kΩC

b
Then the pair (kΩk∞,1 , Ω) of (t, Ω) is feasible. Consequently, the optimality and feasibility of (b
t, Ω)
imply
b ∞,1 ≤ b
bC
b − Ik∞ ≤ 2δ 0 b
t ≤ 2δ 0 kΩk∞,1 .
(54)
kΩk
t ≤ kΩk∞,1 ,
kΩ
Then, on the event E 0 , we obtain

b θb − Ωθ
b + Ωθ
b − βk∞
kβ̄ − βk∞ = kΩ
b ∞,1 kθb − θk∞ + kΩθ
b − βk∞
≤ kΩk
b ∞,1 + kΩCβ
b
≤ δ 0 kΩk
− βk∞
0
b − Ik∞ kβk1
≤ δ kΩk∞,1 + kΩC

bC
b − Ik∞ + kΩ
bC
b − ΩCk
b ∞
≤ δ 0 kΩk∞,1 + kΩ
b ∞,1 kC
b − Ck∞
≤ 3δ 0 kΩk∞,1 + kΩk

(since kβk1 ≤ 1)

≤ 5δ 0 kΩk∞,1 .

The feasibility of βbj implies that kβbj − β̄ j k∞ ≤ µ. By the triangle inequality, we obtain
kβbj − β j k∞ ≤ kβbj − β̄ j k∞ + kβ̄ j − β j k∞ ≤ 2µ,

b = K on the event E gives
since µ = 5δ 0 kΩk∞,1 . Then following from Lemma 12 and using K
bj· − Aj· kq =
kA

K
X
a=1

|βbaj − βaj |q

!1/q



=

X

a∈sj

1/q

|βbaj − βaj |q 

1/q

≤ 2sj µ,

for any 1 ≤ q ≤ ∞. This completes the proof of the last step and of Theorem 4.
To conclude this section we give below the proofs of the intermediary results used in the proof.

37

b = K. Then, from the definition of C
baa ,
Proof of Lemma 11. On the event E, we showed that K
we have
X
1
baa − Caa | ≤ max
b ij | − Caa
max |C
|Σ
1≤a≤K
1≤a≤K |Iba |(|Iba | − 1)
b
i,j∈Ia ,i6=j

E

≤ δ+

1
|Iba |(|Iba | − 1)

X

i,j∈Iba ,i6=j

|Σij | − Caa .

Theorem 3 states that, on the event E, Iba = Ia ∪ La where La = Iba ∩ J1a , for any a ∈ [K]. Therefore,
we consider the following three cases:
(1) For any i, j ∈ Ia and i 6= j, Lemma 7 implies |Σij | − Caa = 0.
(2) For any i ∈ Ia and j ∈ La , the definition of J1a gives

8δ
|Σij | − Caa ≤ 1 − |Aja | (2Caa − ν) ≤ kCk∞ − 4δ.
ν

P
(3) For any i, j ∈ La and i 6= j, since i, j ∈ J1a , we know |Aka | ≥ 1 − 4δ/ν and b6=a |Akb | ≤ 4δ/ν,
for k = i, j. Thus,
X
X
X

|Σij | − Caa ≤ 1 − |Aia ||Aja | Caa + |Aia |
|Ajc ||Cac | + |Aja |
|Aib ||Cab | +
|Aib Ajc |Cbc
c6=a

b6=a

b,c6=a


≤ 1 − |Aia ||Aja | Caa + |Aia |(1 − |Aja |)(Caa − ν) + |Aja |(1 − |Aib |)(Caa − ν)
+ (1 − |Aib |)(1 − |Ajc |)Cb∗ b∗
(for some b∗ 6= a)
"

 #
4δ 2
8δ
16δ 2
≤ 1− 1−
Caa + (Caa − ν) + 2 Cb∗ b∗
ν
ν
ν
≤

16δ
kCk∞ − 8δ,
ν

(by the definition of J1 )

(by (28)).

Therefore, by combining cases (1) - (3), we have
baa − Caa | ≤ δ + |Ia ||La | + |La |(|La | − 1) ·
max |C
1≤a≤K
|Iba |(|Iba | − 1)


16
≤
kCk∞ − 7 δ.
ν



16δkCk∞
− 8δ
ν



where the last inequality comes from that |La | + |Ia | = |Iba |. For the off-diagonal entries, since
bia ) = sign(Aia ), for any i ∈ Ib and a ∈ [K], we have
sign(A
max

1≤a,b≤K,a6=b

bab − Cab | ≤ δ +
|C

we consider the following three cases:
38

1

|Iba ||Ibb |

X

i∈Iba ,j∈Ibb

|Σij | − |Cab | ,

(1) For any i ∈ Ia , j ∈ Ib , we have |Σij | − |Cab | = 0.
(2) For any i ∈ Ia , j ∈ J1b , we have
|Σij | − |Cab | ≤ (1 − |Ajb |)|Cab | +
(3) For any i ∈ J1a , j ∈ J1b , we obtain
Σij

= Aia Ajb Cab + Aia

X

X
c6=b

|Ajc ||Cac | ≤

Ajd Cad +

d6=b

X

8δ
kCk∞ − 4δ.
ν

Aic

c6=a

X

Ajd Ccd .

d∈s(j)

Thus,
≤ (1 − |Aia ||Ajb |)|Cab | + |Aia |(1 − |Ajb |)kCk∞ + (1 − |Aia |)kCk∞




8δ 16δ 2
8δ 16δ 2
≤
− 2
(Caa − ν) +
− 2
kCk∞
ν
ν
ν
ν
16δ
≤
kCk∞ − 8δ.
(by ν < kCk∞ )
ν

|Σij | − |Cab |

Therefore, combining the three cases gives
bab − Cab | ≤ δ + |Ia ||Lb | + |La ||Ib | + 2|La ||Lb | ·
max
|C
1≤a,b≤K,a6=b
2|Iba ||Ibb |


16
kCk∞ − 7 δ.
≤
ν



16δkCk∞
− 8δ
ν



Combining the diagonal and off-diagonal cases yields


16
b
kC − Ck∞ ≤
kCk∞ − 7 δ ≤ 2δ 0
ν

bia ) for any i ∈ I,
b we obtain
We now proceed to bound maxj∈Jb kθbj − θj k∞ . From sign(Aia ) = sign(A
max kθbj − θj k∞ ≤ δ +
j∈Jb

X
1 X
Σij −
Ajb Cab .
a∈[K],j∈Jb |Iba |
max

i∈Iba

P

b∈s(j)

Since for any i ∈ Ia and any j ∈ J, Σij = b∈s(j) Ajb Cab , we focus on the case when i ∈ La . For
any i ∈ La and j ∈ J, (41) yields
X
X
X
X
Ajb Cab − Σij = (1 − Aia )
Ajb Cab −
Aic
Ajb Cbc ,
b∈s(j)

b∈s(j)

c6=a

b∈s(j)

which, by the definition of J1 , implies
X

b∈s(j)

Ajb Cab − Σij

≤ (1 − |Aia |)|Ajd ||Cad | + (1 − |Aia |)|Ajd0 ||Cad0 | (for some d, d0 ∈ [K])
≤

4δ
8δ
(2kCk∞ − ν) ≤
kCk∞ − 4δ.
ν
ν
39

Since we have Jb ⊆ J, we have

|La |
·
max kθbj − θj k∞ ≤ δ + max
a
|Iba |
j∈Jb






8
8
kCk∞ − 4 δ ≤
kCk∞ − 3 δ = δ 0 ,
ν
ν

which concludes the proof of Lemma 11.

Proof of Lemma 12. Let j ∈ Jb be arbitrarily fixed and βbj be the optimal solution of (16) with
µ = 5kΩk∞,1 δ 0 . For simplicity, we remove the super indices. Starting with the following KarushKuhn-Tucker condition:

subject to
we obtain

sign(βba ) + λa sign(βba − β̄a ) = 0,
λa (|βba − β̄a | − µ) = 0,

λa ≥ 0,

for a = {1, . . . , K},





(56)
b
b
b
b
b
0 = sign(βa ) βa − β̄a + λa βa − β̄a = sign(βa ) βa − β̄a + λa µ,

(55)

(56)

(57)

by multiplying both sides of (55) by βba − β̄a . In what follows we prove that if βa = 0, for some a,
then βba = 0. Since this is true when λa = 0 from (55), we only consider when λa 6= 0. Note this
implies |βba − β̄a | = µ from (56). If we assume βba > 0, then (57) gives
β̄a − βba = λa µ.

Since |βba − β̄a | = µ, we further obtain λa = 1 and

β̄a = µ + βba > µ.

(58)

Recall that kβ − β̄k∞ ≤ µ. This implies β̄a ≤ µ + |βa | = µ, which contradicts (58), so βba cannot
be strictly positive. Similarly, βba < 0 cannot hold based on similar arguments. Thus, βba = 0 from
b
which we conclude supp(βbj ) ⊆ supp(β j ) for any j ∈ J.
Proof of Theorem 5. First we construct a set of “hypotheses” of A. Let
M := {v ∈ {0, 1}K : dH (0, v) = s}
where dH (·) denotes the Hamming distance between two binary vectors. Following Lemma A.3 in
Rigollet and Tsybakov [2011] when s ≤ 4K/5, there exists M0 ⊂ M such that, for any w(i) 6=
w(j) ∈ M 0 ,


dH w(i) , w(j) > s/16,
(59)
and

log |M0 | ≥ c0 s log(K/s),

40

(60)

for some constant c0 > 0. We let w(0) = (0, . . . , 0) ∈ RK . Then, we choose
#
"
B
T ∈ Rp×K , for each j = 0, 1, . . . , |M0 |,
A(j) =
η w(j)

where

and




B1


 B2 
(p−1)×K

B= . 
,
∈R
 .. 
BK
η=

r

 
eT
 kT 
ek 
|Ik |×K

Bk = 
, for k ∈ [K],
 ..  ∈ R
.
eTk

c0 c1
8

r

log(1 + K/s)
.
n

(61)

(62)

(63)

2
with some constant c1 = σmin
∧ 1. We use ek to denote the canonical basis of K dimensional space
and 0 to denote the zero vector. Note that, for each Bk , the only non-zero values are at the kth
P
(j) with only one
column. By specifying as above, we choose K
k=1 |Ik | = p − 1 and consider the A
non-pure row. It is easy to verify that A(j) ∈ A for each j = 0, 1, . . . , |M0 | under (31).
We denote by KL(P, Q) the Kullback-Leibler divergence between two probability distributions
P and Q. In order to apply Theorem 2.5 in Tsybakov [2009] to prove (32), for fixed α ∈ (0, 1/8),
we need to check the following three conditions:

(a) KL(PA(i) , PA(0) ) ≤ α log |M0 |, for each i = 1, . . . , |M0 |.
(b) For any 0 ≤ i < j ≤ |M0 |, with some constant c0 > 0,
r


log(K/s)
0 1/q
(i)
(j)
≥cs
.
Lq A , A
n
(c) Lq ( · ) satisfies the triangle inequality.
To show (a), note that X ∼ N (0, ACAT + Γ) from Z ∼ N (0, C = IK ). Invoking Lemma 14
gives


nη 2 s
1
≤
KL PA(i) , PA(0) ≤
log |M 0 |, for any i = 1, . . . , |M 0 |,
(64)
2
16
2σmin

2
by using (60), (63) and c1 = σmin
∧ 1.
To prove (b), for any i = 1, . . . , |M0 |, observe that


Lq A(i) , A(0) = ηkw(i) kq = s1/q η

and, for any i 6= j different from 0,


Lq A(i) , A(j) = ηkw(i) − w(j) kq ≥ (s/16)1/q η ≥ (s1/q η)/16,

by using (59). Combining these two and using the expression of η yield
r


log(K/s)
(i)
(j)
0 1/q
,
Lq A , A
≥cs
n
41

(65)

for 0 ≤ i < j ≤ |M0 |.
e A)
b
Finally, we verify (c) by showing that Lq (·) satisfies the triangle inequality. Consider (A, A,
and observe that
e =
Lq (A, A)

=

≤
=

e ∞,q
min kAP − Ak

P ∈HK

e ∞,q
min kAP − AQk


b ∞,q + kA
b − AQk
e ∞,q
min
kAP − Ak

P,Q∈HK
P,Q∈HK

b ∞,q + min kA
b − AQk
e ∞,q
min kAP − Ak

P ∈HK

Q∈HK

b + Lq (A,
e A).
b
= Lq (A, A)

Therefore, we conclude the proof of (32) by invoking the Theorem 2.5 in Tsybakov [2009].
Lemma 13. (Sherman-Morrison formula) Suppose M ∈ Rn×n is an invertible square matrix and
u, v ∈ Rn are column vectors. Then M + uv T is invertible if and only if 1 + v T M −1 u 6= 0. If
M + uv T is invertible, then its inverse is given by
(M + uv T )−1 = M −1 −

M −1 uv T M −1
.
1 + v T M −1 u

Lemma 14. Assume model (1) and (Z, E) ∼ N (0, diag(IK , Γ)). Let A(0) and A(i) be constructed
as (61) and (62), for any 1 ≤ i ≤ M 0 with M 0 satisfying (60). Let PA(0) and PA(i) be the probability
densities of X parametrized by A(0) and A(i) , respectively. Then we have
KL (PA(i) , PA(0) ) ≤

nη 2 s
2 ,
2σmin

(66)

2
where σmin
= min1≤i≤p Γii .

Proof of Lemma 14. Notationally, we denote by I1 , . . . , IK the set of indices, by Id the d×d diagonal
matrix, by 1|Ik | the |Ik |-dimensional vector and by ΓIk the Ik ×Ik submatrix of Γ. Finally, 0 denotes
the zero-vector and its length might vary line by line.
From the property of Kullback-Leibler divergence, we only need to verify the case when n = 1.
We consider arbitrary A(i) constructed as (61) and (62) for some 0 ≤ i ≤ M 0 . For notational
simplicity, we write A = A(i) = (B T , ξ)T where ξ = ηw(i) ∈ RK . For this A ∈ Rp×K , under
assumptions (a) - (c), observe that
"
#
"
#
T +Γ
BCB
Bξ
Σ
Σ
11
12
B
Σ = ACAT + Γ =
:=
ξT BT
kξk2 + γ
Σ21 Σ22
e 6= A constructed in
by decomposing Γ := diag(ΓB , γ) := diag(ΓI1 , . . . , ΓIK , γ). Similarly, for any A
the same way, we have
"
#
"
#
T +Γ
e
e
BCB
B
ξ
Σ
Σ
11
12
B
e = AC
e A
eT + Γ =
Σ
e 2 + γ := Σ
e 21 Σ
e 22
ξeT B T
kξk
42

Since the Kullback-Leibler divergence between two p-dimensional multivariate normal distributions
N0 := Np (µ0 , Σ0 ) and N1 := Np (µ1 , Σ1 ) is given by




det Σ1
1
> −1
−1
tr Σ1 Σ0 + (µ1 − µ0 ) Σ1 (µ1 − µ0 ) − p + log
,
KL(N0 , N1 ) =
2
det Σ0

we have

 



1
det Σ
−1 e
KL(PAe, PA ) =
tr Σ Σ − p + log
,
e
2
det Σ

(67)

from the fact that X has mean zero. By using the formula of the inverse of a block matrix
"
#−1 "
#
A B
A−1 + A−1 B(D − CA−1 B)−1 CA−1 −A−1 B(D − CA−1 B)−1
=
,
C D
−(D − CA−1 B)−1 CA−1
(D − CA−1 B)−1
for square matrices A and D and non-singular matrices A and D − CA−1 B, we have
#
"
#
"
−1
−1
−1
−1
Ω11 Ω12
Σ−1
−Σ−1
−1
11 + Σ11 Σ12 Σ22·1 Σ21 Σ11
11 Σ12 Σ22·1
:=
,
Σ =
−1
Σ−1
Ω21 Ω22
−Σ−1
22·1
22·1 Σ21 Σ11
with Σ22·1 = Σ22 − Σ21 Σ−1
11 Σ12 . This gives




e 2 + γ) .
e T B T + Ω21 B ξe + Ω22 (kξk
e = tr Ω11 Σ11 + Ω12 (ξ)
tr Σ−1 Σ
|
{z
}
{z
}
|
T
2

T1

We first calculate T1 . Observe that
e T BT
Ω11 Σ11 + Ω12 (ξ)

−1 e T T
−1
−1
= Ip−1 + Σ−1
11 Σ12 Σ22·1 Σ21 − Σ11 Σ12 Σ22·1 (ξ) B

−1
T T
Σ21 = ξ T B T
= Ip−1 + Σ−1
11 Σ12 Σ22·1 ∆ B ,

where ∆ := ξ − ξe ∈ RK . Recall that C = IK . It implies


D1


..
Σ11 = BCB T + ΓB = 
 , with
.
DK

(68)

Dk = 1|Ik | 1T|Ik | + ΓIk .

Thus,

Σ−1
11

 −1
D1

=

..

.
−1
DK




,

with

Dk−1

=

Γ−1
Ik

−

−1
T
Γ−1
Ik 1|Ik | 1|Ik | ΓIk

1 + 1T|Ik | Γ−1
Ik 1|Ik |

, for each k ∈ [K],

P
where the second part uses the Sherman-Morrison formula in Lemma 13. Let αk := i∈Ik 1/σi2 .
Then,
 
eT
!  kT  "
#
−1
−1
T
−1
Γ
1
1
Γ
e


Γ
1
|I
|
|I
|
I
I
k
|I
|
k
k
k
k
 .  = 0 · · · Ik k · · · 0 ∈ R|Ik |×K
Dk−1 Bk = Γ−1
(69)
Ik −
.
1 + αk
1 + αk
.
eTk
43

where only the kth column are non-zero. This further implies
 −1 
 ξ1 −1

D1 B1
1+α1 ΓI1 1|I1 |




..
..
p−1
Σ−1
ξ = 
∈R .
.
.
11 Σ12 = 

(70)

ξK
−1
1+αK ΓIK 1|IK |

−1
DK
BK

On the other hand, note that

h
i
∆T B T = ∆1 1T|I1 | · · · ∆K 1T|IK | ∈ R1×(p−1) .

(71)

Since Σ22·1 is a scalar, combining (68), (70) and (71) yields


−1
T T
T1 = p − 1 + tr Σ−1
11 Σ12 Σ22·1 ∆ B

 ξ1 −1

1+α1 ΓI1 1|I1 |
h
i



..
= p − 1 + Σ−1
 ∆1 1T|I1 | · · · ∆K 1T|IK | 
.
22·1 · tr 
ξK
−1
1+αK ΓIK 1|IK |

= p−1+

Σ−1
22·1

= p−1+

Σ−1
22·1

K

X
ξk ∆k  −1
tr ΓIk 1|Ik | 1T|Ik |
1 + αk

k=1
K
X
k=1

αk
ξk ∆ k .
1 + αk



by using tr

Γ−1
Ik



= αk



Then we proceed to calculate T2 . Similar as (70) and (71), we have

i  ξ
h
ξK
1
−1
−1
T
T
T
−1
−1
T
T
=
∈ R1×(p−1)
=
ξ
1
Γ
·
·
·
1
Γ
Σ21 Σ−1
· · · BK DK
B1 D1
11
1 + α1 |I1 | I1
1 + αK |IK | IK
and




1|I1 | ξe1


B ξe =  ...  ∈ Rp−1 .
1|I | ξeK
K

These two immediately yield

−1 e
−1
e2
T2 = −Σ−1
22·1 Σ21 Σ11 B ξ + Σ22·1 (kξk + γ)

=

Σ−1
22·1

= Σ−1
22·1

K
X
ξk ξek T −1
e +γ−
1 Γ 1
kξk
1 + αk |Ik | Ik |Ik |
k=1
!
K
X
α
k
e 2+γ−
kξk
ξk ξek .
1 + αk
2

!

k=1

Adding up T1 and T2 gives




e = p − 1 + Σ−1
tr Σ−1 Σ
22·1

e 2+γ+
kξk
44

K
X
k=1

!


αk
ξk ∆k − ξek
.
1 + αk

(72)

We proceed to obtain the expression of Σ22·1 . Observe that
 −1
h
i D1
T
..
T 
Σ21 Σ−1
B1T · · · BK

.
11 Σ12 = ξ
= ξT

K
X

BkT Dk−1 Bk

k=1

!

−1
DK




B1
  .. 
 . ξ
BK

ξ

From (69), we further obtain



0
0

 . 
..
 ..  "
.
# 

−1



ΓIk 1|Ik |
 T 
T −1
··· 0 = 
Bk Dk Bk = 1|Ik |  0 · · ·

 . 
1 + αk

 . 

 . 

0


1T
Γ−1 1|Ik |
|I | k
k

1+αk

..

0

which is a K × K diagonal matrix. Thus,

2
Σ22·1 = Σ22 − Σ21 Σ−1
11 Σ12 = kξk + γ −

=

K
X
k=1

.

K
X
k=1





 ∈ RK×K ,





αk
ξ2
1 + αk k

1
ξ 2 + γ.
1 + αk k

(73)

e Recall that the inverse formula of
Finally, it remains to calculate the determinant of Σ and Σ.
a block matrix is
!
A B
det
= det(A) det(D − CA−1 B)
C D
for the invertible matrix A. We thus obtain
e = det Σ11 · Σ
e 22·1 ,
det Σ

det Σ = det Σ11 · Σ22·1 ,

e 22·1 = Σ
e 22 − Σ
e 21 Σ−1 Σ
e
with Σ
11 12 . This gives


det Σ
e 22·1
= log Σ22·1 − log Σ
log
e
det Σ

(74)

Following the same arguments from (73), one can derive
e 22·1 =
Σ

K
X
k=1

Therefore, collecting (72) and (74) yields
2KL PAe, PA



= Σ−1
22·1

e 2+γ+
kξk

K
X
k=1

1 e2
ξ + γ.
1 + αk k



αk
ξk ∆k − ξek
1 + αk
45

(75)

!

e 22·1 ,
− 1 + log Σ22·1 − log Σ

e 22·1 specified in (73) and (75), respectively. Given our specification of A = A(0)
with Σ22·1 and Σ
e = A(i) in (61) and (62), we further have
and A
ξ = 0 =⇒ Σ22·1 = γ,

and

η2s
+ γ,
2
2
1 + Imax /σmin
1 + Imax /σmin

2
= maxk∈[K] |Ik | and σmin
= mini∈[p] σi2 . Consequently, KL PAe, PA can be simplified as
(75)

with Imax

e 2 = η 2 s,
ξe = ηw(i) =⇒ kξk

e 22·1 ≥
Σ

2KL PAe, PA



=

≤

1

e 2+γ =
kξk

η2s
e 22·1
+ log γ − log Σ
γ


η2s
η2s
+γ
+ log γ − log
2
γ
1 + Imax /σmin

≤

η2s
η2s
− 2
2 )γ
γ
η s + (1 + Imax /σmin

=

2
Imax /σmin
+ (η 2 s)/γ
η2s
·
2
γ 1 + Imax /σmin
+ (η 2 s)/γ

≤

η2s
2 ,
σmin

(the concavity of log function)

2 . This concludes the proof of Lemma 14.
by using γ ≥ σmin

A.4

Proofs for the results from Section 4.3

We first prove the three statements of Theorem 6, then present the proofs of Remark 4. Without
loss of generality, we assume that the signed permutation P is identity.
Proof of Theorem 6. We first give the proof for part (a). Then, for ease of the presentation, we
prove part (c) first and then part (b).
b b) ⊆ supp(A b). In addition,
Proof of part (a). Recall that Lemma 12 immediately implies supp(A
J
J
b From the way we construct A
bb, we have supp(A
bb) ⊆
Theorem 3 yields Iba ⊆ Ia ∪ J1a , for any a ∈ K.
I
I
b ⊆ supp(A).
supp(AIb). Therefore, we have proved supp(A)
On the other hand, for any (j, a) ∈ supp(AJ2 ), we know |βaj | > 2µ. This and the fact that
kβbj − β j k∞ ≤ 2µ, immediately gives
|βbaj | ≥ |βaj | − kβbj − β j k∞ > 0,

bJ ).
which implies supp(AJ2 ) ⊆ supp(A
2
bia ) = sign(Aia ) for any (i, a) ∈
b
To show sign(ASb) = sign(ASb), since Lemma 10 guarantees sign(A
b we focus on any fixed (j, a) ∈ Sb and j ∈ J.
b First, we consider the case A
bja = βbaj > 0.
Sb and i ∈ I,
Removing super indices, if βba > 0, (58) gives β̄a > µ. Thus, βa ≥ β̄a − kβ − β̄k∞ > 0 by recalling
bja > 0, (j, a) ∈ Sb and j ∈ J,
b we have Aja > 0.
kβ − β̄k∞ ≤ µ. So far, we have shown that, for any A
b
b
Since the same argument holds for any Aja < 0, the proof of sign(ASb) = sign(ASb) is completed.
46

Proof of part (c). Recall that, for any i ∈ [p] and a ∈ [K],
b a ⇐⇒ A
bia =
i∈G
6 0.

i ∈ Ga ⇐⇒ Aia 6= 0,

We start our proof by rewriting the equivalent expression of TFPP and TFNP:
TFPP =

P

i∈[p],a∈[K] 1{Aia

P

bia 6= 0}
= 0, A

i∈[p],a∈[K] 1{Aia

= 0}

,

TFNP =

P

i∈[p],a∈[K] 1{Aia

P

bia = 0}
6= 0, A

i∈[p],a∈[K] 1{Aia

6= 0}

.

b ⊆ supp(A). Thus,
We first show TFPP = 0. From the result of part (a), we know supp(A)
X
bia 6= 0} = 0,
1{Aia = 0, A
i∈[p],a∈[K]

which implies TFPP = 0.
In order to prove the result of TFNP, observe
X
X
1{Aia 6= 0} = |I| +
si .

(76)

i∈J

i∈[p],a∈[K]

b we partition [p] = I ∪J1 ∪J2 ∪J3 = I ∪(L1 ∪L2 )∪J2 ∪J3
with si = kAi. k0 for each j ∈ J. For given I,
with L1 = Ib∩J1 and L2 = J1 \L1 . Let us consider the set I ∪L1 first. Theorem 3 implies I ∪L1 = Ib
bb, we have
and Iba \ Ia ⊆ J1a . From the way we construct A
I
X

i∈I∪L1 ,a∈[K]

bia = 0} =
1{Aia 6= 0, A

X

i∈L1 ,a∈[K]

bia = 0}.
1{Aia 6= 0, A

Since the definition of J1 implies that, for any j ∈ J1a and a ∈ [K], |Aja | ≥ 1−4δ/ν and |Ajb | ≤ 4δ/ν,
for any b 6= a, this implies
X
X
bia = 0} =
1{Aia 6= 0, A
1{Aib 6= 0} = ti ,
b6=a

a∈[K]

for any i ∈ J1a ∩ L1 and a ∈ [K]. Thus, we have
X
X
bia = 0} =
ti .
1{Aia 6= 0, A

(77)

i∈L1

i∈I∪L1 ,a∈[K]

Next we consider the set L2 . On the event E, for any i ∈ J1a ∩ L2 , we have
bia | ≥ |Aia | − kA
b − Ak∞ ≥ 1 −
|A

bia =
Thus, A
6 0, which implies

X

i∈L2 ,a∈[K]

4δ
− 2µ > 0.
ν

bia = 0} ≤
1{Aia 6= 0, A
47

X

i∈L2

ti .

(78)

bJ ) which yields
Then we consider the set J2 . Part (a) gives supp(AJ2 ) = supp(A
2
X
bia = 0} = 0.
1{Aia 6= 0, A

(79)

i∈J2 ,a∈[K]

bja 6= 0
Finally, we consider the set J3 . By examining the proof of Part (a), it is easy to verify that A
if |Aja | ≥ (2µ) ∨ (4δ/ν), for any j ∈ J3 and a ∈ [K]. Thus,
X
X
bia = 0} ≤
1{Aia 6= 0, A
ti .
(80)
i∈J3

i∈J3 ,a∈[K]

At last, combining (76) - (80) gives
TFNP =

P

i∈[p],a∈[K] 1{Aia

P

bia = 0}
6= 0, A

i∈[p],a∈[K] 1{Aia

6= 0}

≤

P

j∈J1 ∪J3 tj

|I| +

P

j∈J

sj

.

b a ) and GFNP(G
b a ) by the following:
Proof of part (b). Similarly, we can express GFPP(G
ba ) =
GFPP(G

P

i∈[p] 1{Aia

P

bia 6= 0}
= 0, A

i∈[p] 1{Aia

= 0}

ba ) =
GFNP(G

,

P

i∈[p] 1{Aia

P

bia = 0}
6= 0, A

i∈[p] 1{Aia

6= 0}

.

b GFPP(G
b a ) = 0 follows immediately by noting that
For any given a ∈ [K],

ba |
|(Ga )c |
|(Ga )c ∩ G
b a ),
= PK
GFPP(G
0 = TFPP ≥ PK
c|
c|
|(G
)
|(G
)
b
b
b=1
b=1

b a ) = 0 if (Ga )c = 0. To show the expression of GFNP(G
b a ), by the
with the convention GFPP(G
definition of I and Theorem 3, we obtain
X
X
bia = 0} = 0,
1{Aia 6= 0, A
1{Aia 6= 0} = |Ia |.
i∈I

i∈I

The latter immediately implies

X

i∈[p]

1{Aia 6= 0} = |Ia | +

X

sai

i∈J

In addition, following the same arguments in the proof of part (b), we have
X
X
X
bia = 0} = 0,
bia = 0} =
1{Aia 6= 0, A
1{Aia 6= 0, A
i∈J2

i∈J1 ∪J3

tai .

i∈J1 ∪J3 \J1a

Combining these two concludes the proof.

Proofs of Remark 4.
PK
We briefly verify the first claim. It suffices to verify
a=1 |Aja | ≤ 1 which is equivalent with
48

b Recall (16), since β j is feasible, the optimality of βbj immediately gives
kβbj k1 ≤ 1, for any j ∈ J.
kβbj k1 ≤ 1.

To verify the expression of TFNP in the second claim, we assume tj = t and sj = s, for j ∈ J, and
|J1 | + |J3 | = α(|I| + |J3 |). Note that |I| + |J1 | + |J2 | + |J3 | = p implies |J1 | + |J3 | = αp/(1 + α). We
therefore obtain
TFNP ≤
=

t(|J1 | + |J3 |)
t
=
s|J
|+|I|
2
s(|J1 | + |J3 |) + s|J2 | + |I|
s+
· 1+α
α
p

t
by (1 + α)(|I| + |J2 |) = p ,
|I|+s|J
|
s + α1 · |I|+|J22|

as desired.
We verify the third claim. On the event E, when J2 = J, Remark 2 yields Ib = I, Ib = I and
b − Ck∞ ≤ δ and
Jb = J. After careful examination of the proof of Lemma 11, we derive that kC
maxj∈J kθbj − θk∞ ≤ δ, on the event E. Therefore, choosing λ = δ and µ = 3kΩk∞,1 δ proves the
claim, following the proof of Theorems 4 and 6 step by step.
Finally, we verify the fourth claim on the hard-threshold estimator βej for any j ∈ J. For simplicity,
we remove the super indices. Recall that, βe is defined coordinate-wisely by β̄a 1{|β̄a | > µ} with
µ = 5kΩk∞,1 δ 0 .
First, we show kβe − βk∞ ≤ 2µ. For any a ∈ [K] such that |β̄a | ≤ µ, we have
|βea − βa | = |βa | ≤ kβ̄ − βk∞ + |β̄a | ≤ 2µ,

e − Ak∞ ≤ 2µ where
while the same bound is obtained above for the case |β̄a | > µ. This proves kA
e combines A
bb and βej for each j ∈ J.
b To prove the same rate in Theorem 4 for A,
e it suffices to
A
I
j
show that Lemma 12 still holds for βe . Recall that, on the event E, we have kβ̄ − βk∞ ≤ µ. For
any βa = 0, we thus have |β̄a | ≤ kβ̄ − βk∞ ≤ µ, which implies βea = 0. This concludes the proof of
e
Theorem 4 for A.
e and we write (i, a) ∈ Sb if |A
eia | =
To show part (a) of Theorem 6, let Sb denote the support of A
6 0.
b
Let (i, a) ∈ S be arbitrary fixed and consider the following two cases:
b from Theorem 3 and the way we construct A
bb, we have |A
eia | = 1. Thus, |Aia | ≥
- If i ∈ I,
I
eia | − kA
e − Ak∞ ≥ 1 − 2µ > 0.
|A

b then |A
eia | = |βei | = |β̄ i | > µ. Therefore, |Aia | = |β i | ≥ |β̄ i | − kβ̄ i − β i k∞ > 0.
- If i ∈ J,
a
a
a
a

e ⊆ supp(A). To show supp(AJ ) ⊆ supp(A),
e for any (i, a) ∈
Thus, we have proved that supp(A)
2
e
e
supp(AJ2 ), by the definition of J2 , |Aia | > 2µ. Thus, |Aia | ≥ |Aia | − kA − Ak∞ ≥ 0. Therefore,
e
(i, a) ∈ supp(A).
e b) = sign(A b), since Lemma 10 guarantees sign(A
bia ) = sign(Aia ) for any (i, a) ∈
To show sign(A
S
S
b we focus on each (i, a) ∈ Sb and i ∈ J.
b Assuming A
eia = βei > 0, we know βei = β̄ i > µ.
Sb and i ∈ I,
a
a
a
eia > 0 with
Since Aia = βai ≥ β̄ai − kβ̄ i − β i k∞ > 0, we have proved that Aia > 0 for any A
49

b Since the same argument holds for any A
eia < 0, we conclude the proof of
(i, a) ∈ Sb and i ∈ J.
e
sign(ASb) = sign(ASb).
The same conclusion in part (b) and (c) of Theorem 6 holds for GFPP, GFNP, TFPP and
e as it shares the same property in part (a).
TFNP based on the hard-threshold estimator A,

B

Cross-validation illustration

We consider a simple case, when C is diagonal and the signed permutation matrix P is I, to illustrate our cross-validation method.

Example 1. Let C = diag(τ, τ, τ ), I = {1, 2}, {3, 4}, {5, 6} and



1
0
0
 −1
0
0




 0
1
0


 0
1
0

,
A=
0 −1
 0



0 −1
 0


 0.4 0.6 0 
−0.5 0 0.4



∗
τ

0

T
AI CAI = 
0

0
0

τ
∗
0
0
0
0

0
0
∗
τ
0
0

0
0
τ
∗
0
0

0
0
0
0
∗
τ


0
0

0

,
0

τ
∗

where we use ∗ to reflect the fact that our algorithm ignores the diagonal elements. For the true I
bI = AI ,
and I, we have A
b T
b (1) − AI CA
Σ
I
II

For
=
we obtain

F-off

≤

≤


b (1) − ΣII
Σ
II

b (1) − ΣII
Σ
II

b (1) − Σij
max Σ
ij
i6=j



F-off

F-off

b T − ΣII
+ AI CA
I
F-off
p
b
+ |I|(|I| − 1) · kC − Ck∞ .



(2)
b
∨ max Σij − Σij ,
i6=j

1
b (1)
b T
CV (I) = q
 ΣII − AI CAI
|I| |I| − 1

F-off

≤ 2.


Suppose that Ib = {1, 2}, {3, 5}, {4, 6} , so Ib = I, yet Ib =
6 I, we would have

∗ τb1 0 0 0 0
τb1 ∗ 0 0 0 0 


 0 0 ∗ 0 τb 0 


2
bbC
bA
bT = 
A
,
I
Ib
 0 0 0 ∗ 0 τb3 


 0 0 τb2 0 ∗ 0 
0 0 0 τb3 0 ∗




∗
∆τ1 0
∆τ1
∗
0

 0
0
∗
bbC
bA
bT − Σ bb = 
A

I
II
Ib
 0
0
−τ

 0
0
τb2
0
0
0
50

0
0
−τ
∗
0
τb3

0
0
τb2
0
∗
−τ


0
0 

0 

,
τb3 

−τ 
∗

baa from (11) for each a ∈ [K].
b Thus, the
Here ∆τa = τba − τa , using estimates τba defined in lieu of C
cross-validation criterion in (39) would satisfy
1
b b bT
q
 AIbC AIb − ΣIbIb
b |I|
b −1
|I|
s
τ2
4τ 2 + 2b
τ22 + 2b
≥
 3 − 2.
b |I|
b −1
|I|

b ≥
CV (I)

1
b (1)
−q
 ΣIbIb − ΣIbIb
F-off
b |I|
b −1
|I|

F-off

From noting that |b
τa − τ | ≤ , for a = 2, 3, it gives
r
2
2
b ≥ 4τ − 4τ  + 2 − 2 > 2 ≥ CV (I),
CV (I)
15

for τ ≥ 9. We conclude in this example, with Ib = I, incorrectly specifying I will induce a large
b 6= K and Ib =
loss. It is easily verified that this is also the case when Ib = I but K
6 I.
b For instance,
On the other hand, suppose we mistakenly included some non-pure variable in I.

b
suppose we found I = {1, 2}, {3, 4}, {5, 6, 7} . Then we would have




∗
τ
0
0
0 0 0.4τ
∗ τb1 0 0 0 0 0




∗
0
0
0 0 −0.4τ 
 τ
τb1 ∗ 0 0 0 0 0 




 0
 0 0 ∗ τb2 0 0 0 
0
∗
τ
0 0 0.6τ 




bb0 C
bA
bT =  0 0 τb2 ∗ 0 0 0  .
, A
ΣIb0 Ib0 = 
0
0
τ
∗
0
0
0.6τ
0




I
Ib
 0
 0 0 0 0 ∗ τb τb 
0
0
0
∗ τ
0 
3
3







0
0
0
τ ∗
0 
 0
 0 0 0 0 τb3 ∗ τb3 
0.4τ −0.4τ 0.6τ 0.6τ 0 0
∗
0 0 0 0 τb3 τb3 ∗
We thus have



bb0 C
bA
bT − Σ b0 b0
A
I I
I
Ib0

∗

 ∆τ1

 0

=
 0
 0


 0
−0.4τ

∆τ1
∗
0
0
0
0
0.4τ

0
0
∗
∆τ2
0
0
−0.6τ

0
0
∆τ2
∗
0
0
−0.6τ


0
0
−0.4τ

0
0
0.4τ 

0
0
−0.6τ 

0
0
−0.6τ 

∗
∆τ3
τb3 


∆τ3
∗
τb3 
τb3
τb3
∗

and, by similar arguments, for τ ≥ 12, we find
r
4b
τ32 + 4 × 0.36τ 2 + 4 × 0.16τ 2
0
b
CV (I ) ≥
− 2 > 2.
42

Thus, the cross-validation loss in this example will be large even if only one non-pure variable is
mistakenly classified as pure variable. In rare cases, the cross-validation criterion might miss a very
small subset of I but this can be rectified in our later estimation of AJ .

51

C
C.1

Additional Simulation Results
Related work on the estimation of A

As we explained in Section 4.4, the existing procedures for estimating A in (1) are developed
for models satisfying identifiability conditions different than our (i)-(iii). Specifically, Bai and
Li [2012] propose to first optimize, via EM, a quasi-likelihood objective under the identifiability
conditions (a) C = IK and (b) AT Γ−1 A is diagonal. The major advantage of this setting is that the
computationally demanding EM algorithm only needs to determine A and Γ as C = IK is given.
b with the property
The EM algorithm, however, is only guaranteed to find stationary point B
bT B
b is diagonal. In the context of this problem, as the authors note, the EM algorithm
that B
requires a delicate initialization and is computationally demanding, even if only one of K, n and
p is moderately large. Next, the authors propose to link this estimator with an estimator of a
model no longer satisfying (a) and (b) as identifiability conditions, but satisfying instead (1) C is
an arbitrary positive definite matrix; (2) There exists a known set S of K pure variables, with only
one pure variable per latent factor allowed. No further sparsity conditions on A are imposed. To
bB
b T . This yields
estimate A under (2), they suggest to solve for A and C the equation ACAT = B
bS may not be invertible, and
e=B
bB
b −1 of A. However, when K is relatively large, B
the estimator A
S
−1
eS c is not sparse
b = IK , the submatrix A
eS = B
bS B
the estimator may not exist. Finally, although A
S
e
in general. One possibility is to threshold A, but it is unclear how to choose the correct threshold
level, for the following reason. Although the authors establish the asymptotic limit of the MLE of
A under (1) and (2), the estimator of A explained above is not guaranteed to be the MLE in this
model: if it exists, it is a transformation of a stationary point that estimates parameters under the
model specifications (a) and (b), different from (1) and (2). The immediate practical implication
b around A under (1) and (2) is not known, which makes the thresholding
is that the variation of A
b difficult to assess. For all these reasons, we cannot compare numerically our estimation
level of A
procedure with the procedure proposed in Bai and Li [2012], even in the (unrealistic) case when
the pure variable set is known.

C.2

LOVE for non-overlapping cluster estimation

In applications, one may not have prior information on whether the clusters may overlap or not.
Thus, one would prefer a clustering method that works well in both overlapping and non-overlapping
scenarios. In the previous section, we have demonstrated that LOVE outperforms the existing
clustering methods if data are generated from a model that yields variable clusters with overlaps.
In this section, we study the numerical performance of the proposed method under non-overlapping
data generating schemes.
To generate data with non-overlapping clusters, we set the number of variables in each cluster
to be 20. We generate the diagonal elements of C from the uniform distribution in [1, 2] and
use the same method as in Section 5.2 to generate the off-diagonal elements. The variance σj2
of the error Ej is generated from the uniform distribution in [3, 4]. In Table 3, we compare the
sensitivity and specificity of the proposed method with the CORD estimator [Bunea et al., 2016a]
under non-overlapping scenarios, where the sensitivity and specificity are defined in (81). The
52

CORD estimator can be viewed as a benchmark method for variable clustering without overlaps
and is shown to outperform K-means and hierarchical clustering, via an extensive numerical study
presented in Bunea et al. [2016a]. For this reason, we only focus on the comparison between LOVE
and CORD. From Table 3, we see that for small p (i.e., p = 100) the performance of LOVE
is only slightly worse than CORD. As p increases, the specificity of LOVE and that of CORD
remain close to 1, but LOVE yields in fact higher sensitivity than CORD when n = 300. This
confirms that the performance of the proposed method is comparable to the benchmark method
under non-overlapping scenarios. Of course, LOVE is much more flexible as it can detect possible
overlaps.
p

n = 300
LOVE

100
500
1000

n = 500
CORD

LOVE

CORD

SN

SP

SN

SP

SN

SP

SN

SP

0.87
(0.09)
0.86
(0.05)
0.84
(0.05)

0.90
(0.10)
0.98
(0.01)
0.97
(0.02)

0.92
(0.05)
0.82
(0.03)
0.78
(0.03)

0.98
(0.02)
0.98
(0.01)
0.97
(0.01)

0.93
(0.05)
0.87
(0.04)
0.87
(0.04)

0.97
(0.03)
0.99
(0.00)
1.00
(0.00)

0.98
(0.02)
0.94
(0.02)
0.90
(0.02)

1.00
(0.01)
1.00
(0.01)
1.00
(0.01)

Table 3: Sensitivity (SN) and specificity (SP) of the proposed method (LOVE) and CORD under
non-overlapping scenarios. Numbers in parentheses are the simulation standard errors.

C.3

Comparison with other overlapping clustering algorithms

We adopt the same data generating procedure except that we set K = 10 and the negative entries
of A are replaced by their absolute values, since existing overlapping clustering algorithms typically
return an estimator of A with positive entries. We compare the proposed method with the following
overlapping clustering algorithms: fuzzy K-means, and fuzzy K-medoids [Krishnapuram et al.,
2001], the latter being more robust to noise and outliers. We describe the methods briefly in what
follows. Both of them aim to estimate a degree of membership matrix M ∈ Rp×K by minimizing the
ej = (X1j , ..., Xnj ),
average within-cluster L2 or L1 distances [Bezdek, 2013]. Specifically, denote X
n
e = {X
e1 , ..., X
ep }. Let W = {w1 , ..., wK }, where wk ∈ R , be a subset of X
e with K elements.
and X
The fuzzy algorithms aim to find the set W such that J(W ) defined as
J(W ) =

p X
K
X
j=1 k=1

ej , wk ),
Mjk r(X

is minimized. Here, Mjk > 0 can be interpreted as the degree of membership matrix which is a
ej , wk ). Some commonly used expressions of Mjk are shown by Krishnapuram
known function of r(X
53

ej , wk ) is a measure of dissimilarity between X
ej and wk . For instance,
et al. [2001]. In addition, r(X
2
if r(x, θ) = kx−θk2 , this corresponds to the fuzzy K-means. Similarly, the fuzzy K-medoids is given
e is computationally infeasible,
by r(x, θ) = kx − θk1 . Since searching over all possible subsets of X
an approximate algorithm for minimizing J(W ) is proposed by Krishnapuram et al. [2001], we refer
to their original paper for further details.
Their degree of membership matrix M plays the same role as our allocation matrix A, but is
typically non-sparse. In order to construct overlapping clusters based on M one needs to specify a
cut-off value v and assign variable j to cluster k if Mjk > v. Moreover, the number of clusters K
is a required input of the algorithm. In the simulations presented in this section we set K = 10 for
these two methods, which have been implemented by the functions FKM and FKM.med in R.
We compare their performance with our proposed method LOVE. We emphasize that our
method does not require the specification of K and that the tuning parameters are chosen in
a data adaptive fashion, as explained in the previous sections. We follow the pairwise approach of
Wiwie et al. [2015] for this comparison. Recall that G = (G1 , ..., GK ) denotes the true overlapping
b 1 , ..., G
b ) to denote clusters computed from an
clusters. For notational simplicity, we use Gb = (G
K̄
b to be different from K. For
algorithm. Since LOVE estimates the number of clusters, we allow K
any pair 1 ≤ j < k ≤ p, define
n
o
b b for some 1 ≤ a ≤ K and 1 ≤ b ≤ K
b ,
T Pjk = 1 if j, k ∈ Ga and j, k ∈ G
n
o
b b for any 1 ≤ a ≤ K and 1 ≤ b ≤ K
b ,
T Njk = 1 if j, k ∈
/ Ga and j, k ∈
/G
n
o
b b for some 1 ≤ b ≤ K
b ,
F Pjk = 1 if j, k ∈
/ Ga for any 1 ≤ a ≤ K and j, k ∈ G
n
o
b b for any 1 ≤ b ≤ K
b .
F Njk = 1 if j, k ∈ Ga for some 1 ≤ a ≤ K and j, k ∈
/G
and we define
TP =

X

1≤j<k≤p

T Pjk , T N =

X

T Njk , F P =

1≤j<k≤p

X

1≤j<k≤p

F Pjk , F N =

X

F Njk .

1≤j<k≤p

We use sensitivity (SN) and specificity (SP) to evaluate the performance of different methods,
where
TN
TP
SP =
, and SN =
.
(81)
TN + FP
TP + FN
Recall that for the fuzzy methods, variable j belongs to cluster k if the estimated membership
matrix Mjk is beyond a cut-off v, i.e., Mjk > v. We search for the optimal cut-off v in a grid
{0.01, 0.1, ..., 0.3} such that SP + SN is maximized. The corresponding sensitivity and specificity
for LOVE, fuzzy K-means (F-Kmeans) and fuzzy K-medoids (F-Kmed) are shown in Figure 3. To
save space, we only present the results for p = 500 since the other scenarios illustrate the same
patterns. The following findings are observed. First, the F-Kmeans is superior to F-Kmed in most
scenarios in terms of both sensitivity and specificity. Second, LOVE clearly outperforms these
two existing methods and its specificity and sensitivity are very close to 1, which implies that our
method leads to very few false positives and false negatives. The conclusions hold with n from
300 to 1000. Moreover, we reiterate that the true value K = 10 is used as input in the competing
54

●

●

●

1.0

●

●

●

●
●

●

●

0.8

0.8

1.0

methods, whereas it is estimated from the data in LOVE. This illustrates the net advantage of the
proposed method over the existing overlapping clustering methods, for data generated from Model
(1).

●
●

●
●

●
●

●
●

●

0.4

500

700

0.6
0.0

0.0

LOVE
F−Kmeans
F−Kmed

300

0.4

●

0.2

Specificity

0.6

●
●

0.2

Sensitivity

●
●

900

300

sample size

500

700

900

sample size

Figure 3: Plot of specificity and sensitivity for LOVE, fuzzy K-means (F-Kmeans), and fuzzy
K-medoids (F-Kmed) when p = 500.

55

