Provable learning of Noisy-or Networks
Sanjeev Arora∗

Rong Ge†

Tengyu Ma

‡

Andrej Risteski§

arXiv:1612.08795v1 [cs.LG] 28 Dec 2016

December 30, 2016

Abstract

Many machine learning applications use latent variable models to explain structure in
data, whereby visible variables (= coordinates of the given datapoint) are explained as a
probabilistic function of some hidden variables. Finding parameters with the maximum
likelihood is NP-hard even in very simple settings. In recent years, provably efficient
algorithms were nevertheless developed for models with linear structures: topic models,
mixture models, hidden markov models, etc. These algorithms use matrix or tensor
decomposition, and make some reasonable assumptions about the parameters of the
underlying model.
But matrix or tensor decomposition seems of little use when the latent variable model
has nonlinearities. The current paper shows how to make progress: tensor decomposition is applied for learning the single-layer noisy or network, which is a textbook example
of a Bayes net, and used for example in the classic QMR-DT software for diagnosing
which disease(s) a patient may have by observing the symptoms he/she exhibits.
The technical novelty here, which should be useful in other settings in future, is analysis of tensor decomposition in presence of systematic error (i.e., where the noise/error
is correlated with the signal, and doesn’t decrease as number of samples goes to infinity).
This requires rethinking all steps of tensor decomposition methods from the ground up.
For simplicity our analysis is stated assuming that the network parameters were
chosen from a probability distribution but the method seems more generally applicable.

∗

Princeton University, Computer Science Department. arora@cs.princeton.edu
Duke University, Computer Science Department. rongge@cs.duke.edu
‡
Princeton University, Computer Science Department. tengyu@cs.princeton.edu
§
Princeton University, Computer Science Department. risteski@princeton.edu
†

1

Introduction

Unsupervised learning is important and potentially very powerful because of the availability of
the huge amount of unlabeled data — often several orders of magnitudes more than the labeled
data in many domains. Latent variable models, a popular approach in unsupervised learning,
model the latent structures in data: the “structure”corresponds to some hidden variables, which
probabilistically determine the values of the visible coordinates in data. Bayes nets model the
dependency structure of latent and observable variables via a directed graph. Learning parameters
of a latent variable model given data samples is often seen as a canonical definition of unsupervised
learning. Unfortunately, finding parameters with the maximum likelihood is NP-hard even in very
simple settings. However, in practice many of these models can be learnt reasonably well using
algorithms without polynomial runtime guarantees, such as expectation-maximization algorithm,
Markov chain Monte Carlo, and variational inference. Bridging this gap between theory and practice
is an important research goal.
Recently it has become possible to use matrix and tensor decomposition methods to design polynomial-time algorithms to learn some simple latent variable models such as topic models [AGM12, AGH+ 13], sparse coding models [AGMM15, MSS16], mixtures of Gaussians [HK13,
GHK15], hidden Markov models [MR05], etc. These algorithms are guaranteed to work if the
model parameters satisfy some conditions, which are reasonably realistic. In fact, matrix and tensor decomposition are a natural tool to turn to since they appear to be a sweet spot for theory
whereby non-convex NP-hard problems can be solved provably under relatively clean and interpretable assumptions. But the above-mentioned recent results suggest that such methods apply
only to solving latent variable models that are linear: specifically, they need the marginal of the
observed variables conditioned on the hidden variables to depend linearly on the hidden variables.
But many settings seem to call for nonlinearity in the model. For example, Bayes nets in many
domains involve highly nonlinear operations on the latent variables, and could even have multiple
layers. The study of neural networks also runs into nonlinear models such as restricted Boltzmann
machines (RBM) [Smo86, HS06]. Can matrix factorization (or related tensor factorization) ideas
help for learning nonlinear models?
This paper takes a first step by developing methods to apply tensor factorization to learn possibly the simplest nonlinear model, a single-layer noisy-or network. This is a direct graphical model
with hidden variable d ∈ {0, 1}m , and observation node s ∈ {0, 1}n . The hidden variables d1 , . . . , dm
are independent and assumed to have Bernoulli distributions. The conditional distribution Pr[s|d]
is parameterized by a non-negative weight matrix W ∈ Rn×m . We use W i to denote the i-row of
W . Conditioned on d, the observations s1 , . . . , sn are assume to be independent with distribution
Pr [si = 0 | d] =

m
Y

exp(−Wij dj ) = exp(−hW i , di) .

(1.1)

j=1

We see that 1 − exp(−Wji dj ) can be thought of as the probability that dj activates symptom si ,
and si is activated if one of dj ’s activates it — which explains the name of the model, noisy-or. It
follows that the conditional distribution s | d is
Pr[s | d] =

n
Y
i=1

s
1−si
1 − exp(−hW i , di) i exp(−hW i , di)
.

One canonical use of this model is to model the relationship between diseases and symptoms,
as in the classical human-constructed tool for medical diagnosis called Quick Medical Reference
1

(QMR-DT) by (Miller et al.[MPJM82], Shwe et al. [SC91]) This textbook example ([JGJS99]) of
a Bayes net captures relationships between 570 binary disease variables (latent variables) and 4075
observed binary symptom variables, with 45, 470 directed edges, and the Wij ’s are small integers.1
The name “noisy-or ”derives from the fact that Q
the probability that the OR of m independent
binary variables y1 , y2 , . . . , ym is 1 is exactly 1 − j (Pr[yj = 0]). Noisy-or models are implicitly
using this expression; specifically, for the i-th symptom we are considering the OR of m events
where the jth event is “Disease j does not cause symptom i” and its probability is exp(−Wij dj ).
Treating these events as independent leads to expression (1.1).
The parameters of the QMR-DT network were hand-estimated by consulting human experts,
but it is an interesting research problem whether such networks can be created in an automated
way using only samples of patient data (i.e., the s vectors). Previously there were no approaches for
this that work even heuristically at the required problem size (n = 4000). (This learning problem
should not be confused with the simpler problem of infering the latent variables given the visible
ones, which is also hard but has seen more work, including reasonable heuristic methods [JGJS99]).
Halpern et al.[HS13, JHS13] have designed some algorithms for this problem. However, their first
paper [HS13] assumes the graph structure is given. The second paper [JHS13] requires the Bayes
network to be quartet-learnable, which is a strong assumption on the structure of the network.
Finally, the problem of finding a “best-fit” Bayesian network according to popular metrics2 has
been shown to be NP-complete by [Chi96] even when all of the hidden variables are also observed.
Our algorithm and analysis. Our algorithm uses Taylor expansion —on a certain correlation
measure called PMI, whose use in this context is new—to convert the problematic exponential into
an infinite sum, where we can ignore all but the first two terms. This brings the problem into the
realm of tensor decomposition but with several novel twists having to do with systematic error (see
the overview in Section 2). Our recovery algorithm makes several assumptions about W , the matrix
of connection weights, listed in Section 3. We verified some of these on the QMR-DT network, but
the other assumptions are asymptotic in nature. Thus the cleanest description of our algorithm
is in a clean average-case setting. First, we assume all latent variables are iid Bernoulli Ber(ρ)
for some ρ, which should be thought of as small (In the QMR-DT application, ρ is like O(1/m).)
Next we assume that the ground truth W ∈ Rn×m is created by nature by picking its entries in iid
fashion using the following random process:
(
0,
with probability 1 − p
Wij =
f
Wij , with probability p

fij ’s are upper bounded by νu for some constant νu and are identically distributed according
where W
to a distribution D which satisfies that for some constant νl > 0,
h
i
fij2 ) ≤ 1 − νl .
exp(−W
(1.2)
E
fij ∼D
W

fij is bounded away from 0. We will assume that
The condition (1.2) intuitively requires that W
p ≤ 1/3 and νu = O(1), νl = Ω(1). (Again, these are realistic for QMR-DT setting).

1
We thank Randolph Miller and Vanderbilt University for providing the current version of this network for our
research.
2
Researchers resort to these metrics as when the graph structure is unknown, multiple structures may have the
same likelihood, so maximum likelihood is not appropriate.

2

Theorem 1.1 (Informally stated). There exists a polynomial time algorithm (Algorithm 1) that,
given polynomially many samples from the noisy OR network described in the previous paragraph,
e √pm) relative error in ℓ2 -norm in each column.
recovers the weight matrix W with O(ρ

Recall that we mostly thought of the prior of the diseases ρ as being on the order O(1/m). This
√
means that even if p is on the order of 1, our relative error bound equals to O(1/ m) ≪ 1.

2

Preliminaries and overview

We denote by 0 the all-zeroes vector and 1 the all-ones vector. A+ will denote the Moore-Penrose
pseudo-inverse of a matrix A, and for symmetric matrices A, we use A−1/2 as a shorthand for
(A+ )1/2 . The least non-zero singular value of matrix A is denoted σmin (A).
For matrices A, B we define the Kronecker product ⊗ as (A ⊗ B)ijkl = Aij Bkl . A useful identity
is that (A⊗B)·(C ⊗D) = (AC)⊗(BD) whenever the matrix multiplications are defined. Moreover,
Ai will denote the i-th column of matrix A and Ai the i-th row of matrix A.
We write A . B if there exists a universal constant c such that A ≤ cB and we define &
similarly.
The pointwise mutual information of two binary-valued random variables x and y is
E[xy]
P M I2(x, y) , log E[x]
. Note that it is positive iff E[x, y] > E[x] E[y] and thus is used as a
E[y]
measure of correlation in many fields. This concept can be extended in more than one way to a
triple of boolean variables x, y, z and we use
P M I3(x, y, z) , log

E[xy] E[yz] E[zx]
.
[xyz]
E
E[x] E[y] E[z]

(2.1)

(We will sometimes shorten PMI3 and PMI2 to PMI when this causes no confusion.)

2.1

The Algorithm in a nutshell

Our algorithm is given polynomially many samples from the model (each sample describing which
symptoms are or are not present in a particular patient). It starts by computing the following
matrix n × n PMI and and n × n × n tensor PMIT, which tabulate the correlations among all pairs
and triples of symptoms (specifically, the indicator random variable for the symptom being absent):
PMIij , P M I2(1 − si , 1 − sj ).

PMITi,j,k , P M I3(1 − si , 1 − sj , 1 − sj )

(2.2)
(2.3)

The next proposition makes the key observation that the above matrix and tensor are close
to rank m, which we recall is much smaller than n. Here and elsewhere, for a matrix or vector
A we use exp(A) for the matrix or vector obtained by taking the entries-wise exponential. For
convenience, we define F, G ∈ Rn×m as
F , 1 − exp(−W )

G , 1 − exp(−2W ) .

3

(2.4)
(2.5)

Proposition 2.1 (Informally stated). Let Fk , Gk denote the kth columns of the above F, G. Then,


PMI ≈ ρ F F ⊤ + ρGG⊤
PMIT ≈ ρ

m
X
k=1

|

=

{z

m
X

Fk Fk⊤ + ρ2

|

}

m
X
k=1

m
X

Gk G⊤
k

(2.6)

k=1

k=1

Fk ⊗ Fk ⊗ Fk + ρ
:=S

ρ


Gk ⊗ Gk ⊗ Gk .
{z

:=E

(2.7)
(2.8)

}

The proposition is proved later (with precise statement) in Section A by computing the moments
by marginalization and using Taylor expansion to approximate the log of the moments, and ignoring
terms ρ3 and smaller. (Recall that ρ is the probability that a patient has a particular disease,
which should be small, of the order of O(1/n). The dependence of the final error upon ρ appears in
Section 3.) Since the tensor PMIT can be estimated to arbitrary accuracy given enough samples,
the natural idea to recover the model parameters W is to use Tensor Decomposition. This is what
our algorithm does as well, except the following difficulties have to be overcome.
Difficulty 1: Suppose in equation (2.8) we view the first summand S, which is rank m with components Fk ’s as the signal term. In all previous
P polynomial-time algorithms for tensor decomposition,
the tensor is required to have the form m
k=1 Fk ⊗ Fk ⊗ Fk + noise. To make our problem fit this
template we could consider the second summand E as the “noise”, especially since it is multiplied
by ρ ≪ 1 which tends to make E have smaller norm than S. But this is naive and incorrect, since
E is a very structured matrix: it is more appropriate viewed as systematic error. (In particular
this error doesn’t go down in norm as the number of samples goes to infinity.) In order to do tensor
decomposition in presence of such systematic error, we will need both a delicate error analysis and
a very robust tensor decomposition algorithm. These will be outlined in Section 2.3.
Difficulty 2: To get our problem into a form suitable for tensor decomposition requires a whitening
step, which uses the robust estimate of the whitening matrix from the second moment matrix. In
this case, the whitening matrix has to be extracted out of the PMI matrix, which itself suffers from
a systematic error. This also is not handled in previous works, and requires a delicate control of
the error. See Section 2.4 for more discussion.
Difficulty 3: There is another source of inexactness in equation (2.8), namely the approximation
is only true for those entries with distinct indices — for example, the diagonal entry PMIii has
completely different formula from that for PMIij when i 6= j. This will complicate the algorithm,
as described in Subsections 2.3 and 2.4.
The next few Subsections sketch how to overcome these difficulties, and the details appear in
the rest of the paper.

2.2

Recovering matrices in presence of systematic error

In this Section we recall the classical method of approximately recovering a matrix given noisy
estimates of its entries. We discuss how to adapt that method to our setting where the error in
the estimates is systematic and does not go down even with many more samples. The next section
sketches an extension of this method to tensor decomposition with systematic error.
In the classical setting, there is an unknown n × n matrix S of rank m and we are given
S + E where E is an error matrix. The method to recover S is to compute the best rank-m
approximation to S +E. The quality of this approximation was studied by Davis and Kahan [DK70]
and Wedin [Wed72], and many subsequent authors. The quality of the recovery depends upon the
4

ratio ||E||/σm (S), where σm (·) denotes m-th largest singular value and || · || denotes the spectral
norm. To make this familiar lemma fit our setting more exactly, we will phrase the problem as
trying to recover a matrix S given noisy estimate SS ⊤ + E. Now one can only recover S up to
rotation, and the following lemma describes the error in the Davis-Kahan recovery. It also plays a
key role in the error analysis of the usual algorithm for tensor decomposition.
b the subspace of the top m eigenvectors of SS ⊤ and
Lemma 2.2. In the above setting, let K, K
SS ⊤ + E. Let ε be such that kEk ≤ ε · σm (SS ⊤ ). Then IdK − IdKb . ε where Id is the identity
transformation on the subspace in question.
The Lemma thus treats ||E||/σm (SS ⊤ ) as the definition of noise/signal ratio. Before we generalize the definition and the algorithm to handle systematic error it is good to get some intuition,
from looking at (2.6): PMI ≈ ρ(F F T + ρGGT ). Thinking of the first term as signal and the
second as error, let’s check how bad is the noise/signal ratio defined in Davis-Kahan. The “signal”is σm (F F ⊤ ), which is smaller than n since the trace of F F ⊤ is of the order of mn in our
probabilistic model for the weight matrix. The “noise” is the norm of ρGG⊤ , which is large since
the Gk ’s are nonnegative vectors with P
entries of the order of 1, and therefore the quadratic form
h √1n 1, ρGG⊤ √1n 1i can be as large as ρ k hGk , √1n 1i2 ≈ ρmn. Thus the Davis-Kahan noise/signal
ratio is ρm, and so when ρm ≪ 1, it allows recovering the subspace of F with error O(ρm). Note
that this is a vacuous bound since ρ needs to be at least 1/m so that the hidden variable d contains
1 non-zero entry in average. We’ll argue that this error is too pessimistic and we can in fact drive
the estimation error down to close to ρ.
Definition 2.3 (spectral boundedness). Let n ≥ m. Let E ∈ Rn×n be a symmetric matrix and
S ∈ Rn×m . Then, we say E is τ -spectrally bounded by S if
E  τ (SS ⊤ + σm (SS ⊤ ) · Idn )

(2.9)

The smallest such τ is the “error/signal ratio” for this recovery problem.
This definition differs from Davis-Kahan’s because of the τ SS ⊤ term on the right hand side
of (2.9). This allows, for any unit vector x, the quadratic form value xT Ex to be as large as
τ (xT SS ⊤ x + σm (SS ⊤ )). Thus for example the 1 vector no longer causes a large noise/signal ratio
since both quadtratic forms F F ⊤ and GG⊤ have large values on it.
This new error/signal ratio is no larger than the Davis-Kahan ratio, but can potentially be
much smaller. Now we show how to do a better analysis of the Davis-Kahan recovery in terms of
it. The proof of this theorem appears in Section 4.
Theorem 2.4 (matrix perturbation theorem for systematic error). Let n ≥ m. Let S ∈ Rn×m be
of full rank. Suppose positive semidefinite matrix E ∈ Rn×n is ε-spectrally bounded by S ∈ Rn×m
b the subspace of the top m eigenvectors of SS ⊤ and SS ⊤ + E. Then,
for ε ∈ (0, 1). Let K, K
IdK − IdKb . ε .

Finally, we should consider what this new definition of noise/signal ratio achieves. The next
proposition (whose proof appears in Section B) shows that that under the generative model for W
√
sketched earlier, τ = O(log n). Therefore, ρG is Õ(ρ)-bounded by F , and the recovery error of
the subspace of F from F F ⊤ + ρGG⊤ is Õ(ρ) (instead of O(ρm) using Davis-Kahan).
Proposition 2.5. Under the generative model for W , w.h.p, the matrix G = 1 − exp(−2W ) is
τ -spectrally bounded by F = 1 − exp(−W ), with τ = Õ(1).
5

Empirically, we can compute the τ value for the weight matrix W in the QMR-DT
dataset [SC91], which is a textbook application of noisy OR network. For the QMR-DT dataset, τ
is under 6. This implies that the recovery error of the subspace of F guaranteed by Theorem 2.4
is bounded by O(τ ρ) ≈ ρ, whereas the error bound by Davis-Kahan is O(ρm).

2.3

Tensor decomposition with systematic error

Now we extend the insight from the matrix case to tensor recovery under systematic error. In turns
out condition (2.9) is also a good measure of error/signal for the tensor recovery problem of (2.8).
Specifically, if G is τ -bounded by F , then we can recover the components Fk ’s from the PMIT with
√
column-wise error O(ρτ 3/2 m). This requires a non-trivial algorithm (instead of SVD), and the
additional gain is that we can recover Fk ’s individually, instead of only obtaining the subspace with
the PMI matrix.
First we recall the prior state of the art for the error analysis of tensor decomposition with
Davis-Kahan type bounds. The best error bounds involve measuring the magnitude of the noise
matrix Z in a new way. For any n1 × n2 × n3 tensor T , we define the k·k{1}{2,3} norm as
X
xi yjk Tijk .
(2.10)
kT k{1}{2,3} :=
sup
x∈Rn1 ,y∈Rn2 n3
kxk=1,kyk=1

i∈[n1 ]
(j,k)∈[n2 ]×[n3 ]

Note that this norm is in fact the spectral norm of the flattening of the tensor (into a n1 × n2 n3
dimensional matrix). This norm is larger than the injective norm3 , but recently [MSS16] shows that
ε-error in this norm implies O(ε)-error in the recovery guarantees of the components, whereas if one
uses injective norm, the guarantees often pick up an dimension-dependent factor [AGH+ 14]. We
define k·k{2}{1,3} norm similarly. As is customary in tensor decomposition, the theorem is stated
for tensors of a special form, where the components {ui }, {vi }, {wi } are orthonormal families of
vectors. This can be ensured without loss of generality using a procedure called whitening that
uses the 2nd moment matrix.
Theorem 2.6 (Extension of [MSS16, Theorem 10.2]). There is a polynomial-time algorithm (Algorithm 2 later) which has the following guarantee. Suppose tensor T is of the form
T =

r
X
i=1

ui ⊗ vi ⊗ wi + Z

where {ui }, {vi }, {wi } are three collections of orthonormal vectors in Rd , and kZk{1}{2,3} ≤ ε,
kZk{2}{1,3} ≤ ε. Then, it returns {(ũi , ṽi , w̃i )} in polynomial time that is O(ε)-close to {(ui , vi , wi )}
in ℓ2 norm up to permutation. 4
But in our setting the noise tensor has systematic error. An analog of Theorem 2.4 in this
setting is complicated because even the whitening step is nontrivial. Recall also the inexactness
in Proposition 2.1 due to the diagonal terms, which we earlier called Difficulty 3. We address this
difficulty in the algorithm by setting up the problem using a sub-tensor of the PMI tensor. Let
Sa , Sb , Sc be a uniformly random equipartition of the set of indices [n]. Let
ak = Fk,Sa ,
3

bk = Fk,Sb ,

ck = Fk,Sc ,

The injective norm of the tensor T is defined as kT k{1}{2,3} := supx∈Rn1 ,y∈Rn2 z∈Rn3

kxk=1,kyk=1,kzk=1

4

(2.11)
P

i∈[n1 ],j∈[n2 ],k∈[n3 ]

xi yj zk Tijk .

Precisely, here we meant that there exists a permutation π such that for every i, max{kũπ(i) − ui k, kṽπ(i) −
vi k, kw̃π(i) − wi k} ≤ O(ε)

6

where Fk,S denotes the restriction of vector Fk to subset S. Moreover, let
γk = Gk,Sa ,

δk = Gk,Sb ,

θk = Fk,Sc .

(2.12)

Then, since the sub-tensor PMITSa ,Sb ,Sc only contains entries with distinct indices, we can use
Taylor expansion (see Lemma A.1) to obtain that
X
X
γk ⊗ δk ⊗ θk + higher order terms .
PMITSa ,Sb ,Sc = ρ
ak ⊗ bk ⊗ ck + ρ2
k∈[m]

k∈[m]

Here the second summand on the RHS corresponds to the second order term in the Taylor expansion.
It turns out that the higher order terms are multiplied by ρ3 and thus have negligible Frobenius
norm, and therefore discussion below will focus on the first two summands.
For simplicity, let T = PMITSa ,Sb ,Sc . Our goal is to recover the components ak , bk , ck from the
approximate low-rank tensor T .
The first step is to whiten the components ak ’s, bk ’s and ck ’s. Recall that ak = Fk,Sa is a
non-negative vector. This implies the matrix A = [a1 , . . . , am ] must have a significant contribution
in the direction of the vector 1, and thus is far away from being well-conditioned. For the purpose
of this section, we assume for simplicity that we can access the covariance matrix defined by the
vector ak ’s,
X
(2.13)
Q̄a := AA⊤ =
ak a⊤
k .
k∈[m]

Similarly we assume the access of Q̄b and Q̄c which are defined analogously. In Section 2.4 we
discuss how to obtain approximately these three matrices.
Then,
we can compute the whitened tensor by applying transformation
1/2 , (Q̄+ )1/2 , (Q̄+ )1/2 along the three modes of the tensor T ,
)
(Q̄+
c
a
b
1/2
1/2
1/2
⊗ (Q̄+
·T =ρ
(Q̄+
⊗ (Q̄+
c )
a)
b )

k∈[m]

+ ρ2
|

X

1/2
1/2
1/2
bk ⊗ (Q̄+
ck
(Q̄+
ak ⊗ (Q̄+
c )
a)
b )

X

1/2
1/2
1/2
(Q̄+
γk ⊗ (Q̄+
δk ⊗ (Q̄+
θk +negligible terms
a)
c )
b )

k∈[m]

{z

}

:=Z

1/2 a ’s are orthonormal
Now the first summand is a low rank orthogonal tensor, since (Q̄+
k
a)
vectors. However, the term Z is a systematic error and we use the following Lemma to control its
k·k{1}{2,3} norm.

Lemma 2.7. Let n ≥ m and A, B, C ∈ Rn×m be full rank matrices and let Γ, ∆, Θ ∈ Rd×ℓ . Let
γi , δi , θi be the i-th column of Γ, ∆, Θ, respectively. Let Q̄a = AA⊤ , Q̄b = BB ⊤ , Q̄c = CC ⊤ . Suppose
ΓΓ⊤ (and ∆∆⊤ , ΘΘ⊤ ) is τ -spectrally bounded by A (and B, C respectively), then,
X

1/2
1/2
1/2
(Q̄+
γi ⊗ (Q̄+
δi ⊗ (Q̄+
θi
c )
a)
b )

i∈[ℓ]

≤ (2τ )3/2 .
{1}{2,3}

Lemma 2.7 shows that to give an upper bound on the k·k{1}{2,3} norm of the error tensor Z,
it suffices to show that the square of the components of the error, namely, ΓΓ⊤ , ∆∆⊤ , ΘΘ⊤ are
7

τ -spectrally bounded by the components of the signal A, B, C respectively. This will imply that
kZk{1}{2,3} ≤ (2τ )3/2 ρ2 .
Recall that A and Γ are two sub-matrices of F and G. We have shown that GG⊤ is τ -spectrally
bounded by F in Proposition 2.5. It follows straightforwardly that the random sub-matrices also
have the same property.
Proposition 2.8. In the setting of this section, under the generative model for W , w.h.p, we have
that ΓΓ⊤ is τ -spectrally bounded by A with τ = O(log n). The same is true for the other two
modes.
Using Proposition 2.8 and Lemma 2.7, we have that
kZk{1}{2,3} . ρ2 log3/2 (n) .
1/2 ⊗ (Q̄+ )1/2 ⊗ (Q̄+ )1/2 · T , we can recover the comThen using Theorem 2.6 on the tensor (Q̄+
c
a)
b
+ 1/2
+
1/2
+
1/2
ponents (Q̄a ) ak ’s, (Q̄b ) bk ’s, and (Q̄c ) ck ’s. This will lead us to recover ak ,bk and ck , and
finally to recover the weight matrix W .

2.4

Robust whitening

In the previous subsection, we assumed the access to Q̄a , Q̄b , Q̄c (defined in (2.13)) which turns out
to be highly non-trivial. A priori, using equation (2.6), noting that A = [F1,Sa , . . . , Fm,Sa ], we have
PMISa ,Sa /ρ ≈ Q̄a + error .
However, this approximation can be arbitrarily bad for the diagonal entries of PMI since equation (2.6) only works for entries with distinct indices. (Recall that this is why we divided the
indices set into Sa , Sb , Sc and studied the asymmetric tensor in the previous subsection). Moreover,
the diagonal of the matrix Q̄a contributes to its spectrum significantly and therefore we cannot get
meaningful bounds (in spectral norm) by ignoring the diagonal entries.
This issue turns out to arise in most of the previous tensor papers and the solution was to
compute AA⊤ by using the asymmetric moments AB ⊤ , BC ⊤ , CA⊤ ,
AA⊤ = (AB ⊤ )(CB ⊤ )+ (CA⊤ ) .
Typically AB ⊤ , BC ⊤ , CA⊤ can be estimated with arbitrarily small error (as number of samples
go to infinity) and therefore the equation above leads to accurate estimate to AA⊤ . However,
in our case the errors in the estimate PMISa ,Sb ≈ AB ⊤ , PMISb ,Sc ≈ BC ⊤ , PMISc ,Sa ≈ CA⊤ are
systematic. Therefore, we need to use a more delicate analysis to control how the error accumulates
in the estimate,
Q̄a ≈ PMISa ,Sb · PMI−1
Sb ,Sc · PMISc ,Sa .
Here again, to get an accurate bound, we need to understand how the error in PMISa ,Sb −
AB ⊤ behaves relatively compared with AB ⊤ in a direction-by-direction basis. We generalized
Definition 2.3 to capture the asymmetric spectral boundedness of the error by the signal.
Definition 2.9 (Asymmetric spectral boundedness). Let n ≥ m and B, C ∈ Rn×m . We say a
matrix E ∈ Rn×n is ε-spectrally bounded by (B, C) if E can be written as:
⊤
E = B∆1 C ⊤ + B∆⊤
2 + ∆3 C + ∆4 .

(2.14)

Here ∆1 ∈ Rm×m , ∆2 , ∆3 ∈ Rn×m and ∆4 ∈ Rn×n are matrices whose spectral norms are bounded
by: k∆1 k ≤ ε, k∆2 k ≤ εσmin (C), k∆3 k ≤ εσmin (B) and k∆4 k ≤ εσmin (B)σmin (C).
8

Let K be the column subspace of B and H be the column subspace of C. Then we have
∆1 = B + E(C ⊤ )+ , ∆2 = B + EIdH ⊥ , ∆3 = IdK ⊥ E(C ⊤ )+ , ∆4 = IdK ⊥ EIdH ⊥ . Intuitively, they
measure the relative relationship between E and B, C in different subspaces. For example, ∆1 is
the relative perturbation in the column subspace of K and row subspace of H. When B = C,
this is equivalent to the definition in the symmetric setting (this will be clearer in the proof of
Theorem 2.4).
Theorem 2.10 (Robust whitening theorem). Let n ≥ m and A, B, C ∈ Rn×m .
Σab , Σbc , Σca ∈ Rn×n are of the form,

Suppose

Σab = AB ⊤ + Eab , Σbc = BC ⊤ + Ebc , and Σca = CA⊤ + Eca .
where Eab , Ebc , Eca are ε-spectrally bounded by (A, B), (B, C), (C, A) respectively. Then, the matrix
matrix
+
Qa = Σab [Σ⊤
bc ]m Σca
+
⊤
is a good approximation of AA⊤ in the sense that Qa = Σab [Σ⊤
bc ]m Σca − AA is O(ε)-spectrally
bounded by A. Here [Σ]m denotes the best rank-m approximation of Σ.

The theorem is non-trivial even if the we have an absolute error assumption, that is, even
if kEbc k ≤ τ σmin (B)σmin (C), which is stronger condition than Ebc is τ -spectrally bounded by
⊤ +
⊤
(B, C). Suppose we establish bounds on kΣab − AB ⊤ k, kΣ+⊤
bc − (BC ) k and kΣab − AB k
+
⊤
individually, and then putting them together in the obvious way to control the error Σab [Σbc ]m Σca −
AB ⊤ (BC ⊤ )+ CA⊤ . Then the error will be too large for us. This is because standard matrix

⊤ −1
⊤ −1 2 .
perturbation theory gives that kΣ−⊤
bc − (BC ) k can be bounded by O kEbc kk(BC ) k
ε/[σmin (B)σmin (C)], which is tight. Then we multiply the error with the norm of the rest of the
(B)σmax (C)
. That is, we will loss a condition number of
two terms, the error will be roughly ε · σσmax
min (B)σmin (C)
B, C, which can be dimension dependent for our case.
+
The fix to this problem is to avoid bounding each term in Σab [Σ⊤
bc ]m Σca individually. To do
this, we will take the cancellation of these terms into account. Technically, we re-decompose the
+
+
+
⊤ +
product Σab [Σ⊤
bc ]m Σca into a new product of three matrices (Σab B )(B[Σbc ]m C)(C Σca ), and then
bound the error in each of these terms instead. See Section C for details.
1/2 a ’s are indeed approximately
As a corollary, we conclude that the whitened vectors (Q+
i
a)
orthonormal.
1/2 A contains approximately
Corollary 2.11. In the setting of Theorem 2.10, we have that (Q+
a)
orthonormal vectors as columns, in the sense that
1/2
1/2
k(Q+
AA⊤ (Q+
− Idk . ε .
a)
a)

Therefore we have found an approximate whitening matrix for A even though we do not have
access to the diagonal entries.

3

Main Algorithms and Results

As sketched in Section 2, our main algorithm (Algorithm 1) uses tensor decomposition on the PMI
tensor. In this section, we describe the different steps and how the fit together. Subsequently, all
steps will be analyzed in separate sections.

9

Algorithm 1 Learning Noisy-Or Networks via Decomposing PMI Tensor
Inputs: N samples generated from a noisy-or network, disease prior ρ
c.
Outputs: Estimate of weight matrix W

d PMIT
\ using equation (E.1).
1. Compute the empirical PMI matrix and tensor PMI,
2. Choose a random equipartition Sa , Sb , Sc of [n].

\ via Algorithm 4 for the partitioning
3. Obtain approximate whitening matrices for PMIT
Sa , Sb , Sc
4. Run robust tensor-decomposition Algorithm 3 to obtain vectors âi , b̂i , ĉi , i ∈ [m]
1/3 â , 1 − ( 1−ρ )1/3 b̂ , 1 − ( 1−ρ )1/3 ĉ .
5. Let Yi be the concatenation of the three vectors 1 − ( 1−ρ
i
i
i
ρ )
ρ
ρ

(Recall that âi , b̂i , ĉi are of dimension n/3 each.)

c , where
6. Return W

ci,j
W

(
− log((Yi )j ), if (Yi )j > exp(−νu ) :
=
exp(−νu ),
otherwise

Theorem 3.1 (Main theorem, random weight matrix). Suppose the true W is generated from
the random model in Section 1 with ρpm ≤ c for some sufficiently small constant c. Then given
c in polynomial
N = poly(n, 1/p/, 1/ρ) number of examples, Algorithm 1 returns a weight matrix W
time that satisfies
ci − Wi k2 ≤ O(η
e √pn) ,
∀i ∈ [m], kW

√
where η = Õ mpρ .
√
Note that the column ℓ2 norm of Wi is on the order of pn, and thus η can be thought of as the
relative error in ℓ2 norm. Note also that Pr[si = 0] = 1 − Pr[si = 1] ≈ 1 − pmρ, so ρpm = o(1) is
necessary purely for sample complexity reasons. Finally, we can also state a result with a slightly
weaker guarantee, but with only deterministic assumptions on the weight matrix W . Recall that
F = 1 − exp(−W ) and G = 1 − exp(−2W ). We will also define third and fourth-order terms
H = 1 − exp(−3W ), L = 1 − exp(−4W ).
We also define the incoherence of a matrix F . Roughly speaking, it says that the left singular
vectors of F don’t correlate with any of the natural basis vector much more than the average.
Definition 3.2 (Incoherence:). Let Fp∈ Rn×m have singular value decomposition F = U ΣV ⊤ . We
say F is µ-incoherent if maxi kUi k ≤ µm/n. where Ui is the i-th row of U .

We assume the weight matrix W satisfies the following deterministic assumptions,
1. GG⊤ , HH ⊤ , LL⊤ is τ -spectrally
p bounded by F for τ ≥ 1.
e n/m).
2. F is µ-incoherent with µ ≤ O(
3. If maxi kFi k0 ≤ pn, with high probability over the choice of a subset Sa , |Sa | = n/3, σmin (FSa ) &
√
np and ρpm ≤ c for some sufficiently small constant c.
Theorem 3.3 (Main theorem, deterministic weight matrix). Suppose the matrix W satisfies the
c in polynomial
conditions 1-3 above. Given polynomial number of samples, Algorithm 1 returns W
time, s.t.
ci − Wi k2 ≤ O(η
e √np) .
∀i ∈ [m], kW
10

for η =

√

mpρτ 3/2

√
√
Since the ℓ2 norm of Wi is on the order of np, the relative error in ℓ2 -norm is as most mρτ 3/2 ,
which mirrors the randomized case above.
The proofs of Theorems uses the overall strategy of Section 2, and is deferred to Section D. We
give a high level outline that demonstrates how the proofs depends on the machinery built in the
subsequent sections.
Both Theorem 3.1 and Theorem 3.3 are similarly proved – the only technical difference being
how the third and higher order terms are bounded. (Because of generative model assumption, for
Theorem 3.1 we can get a more precise control on them.) Hence, we will not distinguish between
them in the coming overview.
Overall, we will follow the approach outlined in Section 2. Let us step through Algorithm 1 line
by line:
1. The overall goal will be to recover the leading terms of the PMI tensor. Of course, we get
samples only, so can merely get an empirical version of it. In Section E, we show that the
simple plug-in estimator does the job – and does so with polynomially many samples.
2. Recall Difficulty 3 from Section 2 : the PMI tensor and matrix expression is only accurate
on the off-diagonal entries. In order to address this, in Section 2.3 we passed to a sub-tensor
of the original tensor by partitioning the symptoms into three disjoint sets, and considering
the induced tensor by this partition.
3. In order to apply the robust tensor decomposition algorithm from Section 5, we need to first
calculate whitening matrices. This is necessarily complicated by the fact that the diagonals
of the PMI matrix are not accurate, as discussed in Section 2.4. Section C gives guarantees
on the procedure for calculating the whitening matrices.
4. This is main component of the algorithm: the robust tensor decomposition machinery. In
Section 5, the conditions and guarantees for the success of the algorithm are formalized.
There, we deal with the difficulties layed out in Section 2.2 : namely that we have a substantial
systematic error that we need to handle. (Both due to higher-order terms, and due to the
missing diagonal entries)
5. This step, along with Step 6, is a post-processing step – which allows us to recover the weight
matrix W after we have recovered the leading terms of the PMI tensor.
We also give a short quantitative sense of the guarantee of the algorithm. (The reader can find
the full proof in Section D.)
To get quantitative bounds, we will first need a handle on spectral properties of the random
model: these are located in Section B. As we mentioned above, the main driver of the algorithm is
step 4, which uses our robust tensor decomposition machinery in Section 5. To apply the machinery,
we first need to show that the second (and higher) order terms of the PMI tensor are spectrally
bounded. This is done by applying Proposition B.4, which roughly shows the higher-order terms are
O(ρ log n)-spectrally bounded by ρF F ⊤ . The whitening matrices are calculated using machinery
in Section C. We can apply these tools since the random model gives rise to a O(1)-incoherent F
matrix as shown in Lemma B.1.
To get a final sense of what the guarantee is, the l2 error which step 4 gives, via Theorem 5.4
√
roughly behaves like σmax τ 3/2 , where σmax is the spectral norm of the whitening matrices and
τ is the spectral boundedness parameter. But, by Lemma D.1 σmax is approximately the spectral
11

norm of ρF F ⊤ – which on the other hand by Lemma B.1 is on the order of mnp2 ρ. Plugging in
these values, we get the theorem statement.

4

Finding the Subspace under Heavy Perturbations

In this section, we show even if we perturb a matrix SS ⊤ with an error whose spectral norm might
be much larger than σmin (SS ⊤ ), as long as E is spectrally bounded the top singular subspace of
S is still preserved. We defer the proof of the asymmetric case (Theorem 2.10) to Section C. We
note that such type of perturbation bounds, often called relatively perturbation bounds, have been
studied in [Ips98, Li98a, Li98b, Li97]. The results in these papers either require the that signal
matrix is full rank, or the perturbation matrix has strong structure. We believe our results are new
and the way that we phrase the bound makes the application to our problem convenient. We recall
Theorem 2.4, which was originally stated in Section 2.
Theorem 2.4 (matrix perturbation theorem for systematic error). Let n ≥ m. Let S ∈ Rn×m be
of full rank. Suppose positive semidefinite matrix E ∈ Rn×n is ε-spectrally bounded by S ∈ Rn×m
b the subspace of the top m eigenvectors of SS ⊤ and SS ⊤ + E. Then,
for ε ∈ (0, 1). Let K, K
IdK − IdKb . ε .

Proof. We can assume ε ≤ 1/10 since otherwise the statement is true (with a hidden constant 10).
Since E is a positive semidefinite matrix, we write E = RR⊤ where R = E 1/2 . Since A has full
column rank, we can write R = AS + B where S ∈ Rm×n and the columns of B are in the subspace
K ⊥ . (Specifically, we can choose S = A+ R and B = R − AA+ R = IdK ⊥ B.) By the definition of
spectral boundedness, we have


⊤
⊤
⊤
⊤
BB = IdK ⊥ RR IdK ⊥  IdK ⊥ ε AA + σm (AA )Idn IdK ⊥ .
= εσm (AA⊤ )IdK ⊥ .

Therefore, we have that kBk2 ≤ εσmin (AA⊤ ). Moreover, we also have
IdK RR⊤ IdK  εAA⊤ + εσmin IdK ,
It follows that
ASS ⊤ A⊤ ≤ 2εAA⊤ .
which implies
Let P = Idm + SS ⊤

1/2

kSS ⊤ k ≤ ε.
. Then we write AA⊤ + E as,

AA⊤ + E = AA⊤ + RR⊤ = AA⊤ + (AS + B)(AS + B)⊤
= A(Id + SS ⊤ )A⊤ + ASB ⊤ + BS ⊤ A⊤ + BB ⊤
= (AP + BS ⊤ P −1 )(AP + BS ⊤ P −1 )⊤ + BB ⊤ − BS ⊤ P −2 SB ⊤

(4.1)

b = (AP + BS ⊤ P −1 ). Let K ′ be the column span of A.
b We first prove that K
b is close to K ′ .
Let A
Note that
BB ⊤ − BS ⊤ P −2 SB ⊤ . kBk2 + kBk2 S ⊤ P −2 S . kBk2
. εσmin (AA⊤ ) .
12

(since P = Id + SS ⊤  SS ⊤ )


bA
b⊤ ) = σmin (A)
b 2 = σmin (AP ) − kBS ⊤ P −1 k 2 ≥ (1 − O(ε))σmin (A)2 .
Moreover, we have σmin (A
Therefore, using Wedin’s Theorem (Lemma F.2) on equation (4.1), we have that
kIdKb − IdK ′ k . ε .

Next we show K ′ and K are also close. We have
p
b − AP k ≤ kBS ⊤ P −1 k ≤ ε σmin (A)2
kA

(4.2)

(since kSk .

√

ε, kBk .

√

ε)

b is close to the
Therefore, by Wedin’s Theorem, K ′ , as the span of top m left singular vectors of A,
span of the top left singular vector of AP , namely, K
kIdK − IdK ′ k . ε .

(4.3)

Therefore using equation (4.2) and (4.3) and triangle inequality, we complete the proof.

5

Robust Tensor Decomposition with Systematic Error

In this section we discuss how to robustly find the tensor decomposition even in presence of systematic error. We first illustrate the main techniques in an easier setting of orthogonal tensor
decomposition (Section 5.1), then we describe how it can be generalized to the general setting that
we require for our algorithm (Section 5.2).

5.1

Warm-up: Approximate Orthogonal Tensor Decomposition

We start with decomposing an orthogonal tensor with systematic error. The algorithm we use here
is a slightly more general version of an algorithm in [MSS16].
Algorithm 2 Robust orthogonal tensor decomposition
Inputs: Tensor T ∈ Rd×d×d , number δ, ε ∈ (0, 1).
Outputs: Set S = {(ãi , b̃i , c̃i )}
1. S = ∅
2. For s = 1 to O(d1+δ log d)
3.
4.
5.


Draw g ∼ N (0, Idn ), and compute M = Idn ⊗ Idn ⊗ g⊤ · T .5

Compute the top left and right singular vectors u, v ∈ Rd of M . Let z = (u⊤ ⊗v ⊤ ⊗Idn )·T .

If (u⊤ ⊗ v ⊤ ⊗ z ⊤ ) · T ≥ 1 − ζ, where ζ = O(ε), and u is 1/2-far away from any of ui ’s
with (ui , vi , wi ) ∈ S, then add (u, v, w) to S.

6. Return S

Theorem 5.1 (Stronger version of Theorem 2.6). Suppose {ui }, {vi }, {wi } are three collection
ε-approximate orthonormal vectors. Suppose tensor T is of the form
T =

r
X
i=1

5

ui ⊗ vi ⊗ wi + Z

Recall that product of two tensor (A ⊗ B ⊗ C) · (E ⊗ D ⊗ F ) = AE ⊗ BD ⊗ CF

13

with kZk{2}{1,3} ≤ τ and kZk{1}{2,3} ≤ τ . Then, with probability at least 0.9, Algorithm 2 returns
S = {(ũi , ṽi , w̃i )} which is guaranteed to be O((τ + ε)/δ)-close to {(ui , vi , wi )} in ℓ2 -norm up to
permutation.
Proof Sketch of Theorem 5.1. The Theorem is a direct extension of [MSS16, Theorem 10.2] to
asymmetric and approximate orthogonal case. We only provide a proof sketch here. We start by
writing
m


X
⊤
hg, wi iui vi⊤ + (Idn ⊗ Idn ⊗ g⊤ ) · Z
·T =
M = Id ⊗ Id ⊗ g
|
{z
}
i=1
:=Mg
{z
}
|

(5.1)

:=Ms

Since kZk{2}{1,3} ≤ τ and kZk{1}{2,3} ≤ τ , [MSS16, Theorem 6.5] implies that with probability at
least 1 − d2 over the choice of g,
p
(Idn ⊗ Idn ⊗ g⊤ ) · Z ≤ 2 log d · τ

√
Let t = 2 log d. We have that with probability 1/(d1+δ logO(1) d), hg, w1 i ≥ (1 + δ/3)t and
hg, wj i ≤ t for every j 6= 1. We condition on these events. Let ūi be a set of orthonormal vectors
such that Eu = [u1 , . . . , um ] − [ū1 , . . . , ūm ] satisfies kEu k ≤ ε (we can take ūi ’s to be the whitening
of uP
i ’s). Similarly define v̄i ’s. Then we have that the
Pterm (defined in equation (5.1)) can be written
′
′
as i hg, w1 iūi v̄i + E where kEk . ε. Let M̄S = i hg, w1 iūi v̄i . Then M̄S has top singular value
hg, w1 i ≥ (1 + δ/3)t, and second singular value at most t. Moreover, the term Mg + E ′ has spectral
norm bounded by O(τ + ε). Thus by Wedin’s Theorem (Lemma F.2), the top left and right singular
vectors u, v of MS + Mg = M̄S + Mg + E ′ are O((τ + ε)/δ)-close to ū1 and v̄1 respectively. They
are also O((τ + ε)/δ)-close to u1 , v1 since u1 is close to ū1 . Moreover, we have (u⊤ ⊗ v ⊤ ⊗ Id) · T is
O(τ /δ)-close to w1 .
Therefore, with probability 1/(d1+δ logO(1) d), each round of the for loop in Algorithm 2 will find
u1 , v1 , w1 . Line 5 is used to verify if the resulting vectors are indeed good using the injective norm
as a test. It can be shown that if the test is passed then (u, v, z) is close to one of the component.
Therefore, after d1+δ logO(1) d iterations, with high probability, we can find all of the components.

5.2

General tensor decomposition

In many previous works, general tensor decomposition is reduced to orthogonal tensor decomposition via a whitening procedure. However, here in our setting we cannot estimate the exact whitening
matrix because of the systematic error. Therefore we need a more robust version of approximate
whitening matrix, which we define below:
Definition 5.2. Let r ≤ d. A collection of r vectors {a1 , . . . , ar } is ε-approximately orthonormal if
the matrix A with ai as columns satisfies
A⊤ A − Id ≤ ε

(5.2)

Definition 5.3. Let d ≥ r and A = [a1 , . . . , ar ] ∈ Rd×r . A PSD matrix Q ∈ Rd×d is an εapproximate whitening matrix for A if (Q+ )1/2 A is ε-approximately orthonormal .
With this in mind, we can state the guarantee on the tensor decomposition algorithm (Algorithm 3).
14

Algorithm 3 Tensor decomposition with systematic error
Inputs: Tensor T ∈ Rn1 ×n2 ×n3 and ε-approximate whitening matrices Qa , Qb , Qc ∈ Rd×d .
Outputs: {âi , bˆi , ĉi }i∈[r]
1/2 ⊗ (Q+ )1/2 ⊗ (Q+ )1/2 · T
1. Compute T̃ = (Q+
c
a)
b

2. Run orthogonal tensor decomposition (Algorithm 2) with input T̃ , and obtain {ăi , b̆i , c̆i }
1/2

1/2

1/2

3. Return: {Qa ăi , Qb b̆i , Qc c̆i }

Theorem 5.4. Let d ≥ r, and A, B, C ∈ Rd×r be full rank matrices. Let Γ, ∆, Θ ∈ Rd×ℓ . Let
ai , bi , ci , γi , δi , θi be the columns of A, B, C, Γ, ∆, Θ respectively. Suppose tensor T is of the form
T =

r
X
i=1

ai ⊗ bi ⊗ ci +

ℓ
X
i=1

γi ⊗ δi ⊗ θi + E

(5.3)

Suppose matrices Qa ∈ Rd×d , Qb ∈ Rd×d , Qc ∈ Rd×d are ε-approximate whitening matrices for
1/2
1/2
1/2
A, B, C, and suppose Γ, ∆, Θ are τ -spectrally bounded by Qa , Qb , Qc , respectively. Then,
Algorithm 3 returns âi , b̂i , ĉi that are O(η)-close to ai , bi , ci in Õ(d4+δ ) time with


η . max(kQa k , kQb k , kQc k)1/2 · τ 3/2 + σ −3/2 kEk{1,2}{3} + ε · 1/δ

where σ = min(σmin (Qa ), σmin (Qb ), σmin (Qc )).

Note that in our model, the matrix E has very small spectral norm as it is the third order term
in ρ (and ρ = O(1/n)). The spectral boundedness of Γ, ∆, Θ are discussed in Section B. Therefore
we can expect the RHS to be small.
In order to prove this theorem, we show after we apply whitening operation using the approximate whitening matrices, the tensor is still close to an orthogonal tensor. To do that, we need the
following lemma which is a useful technical consequence of the condition (2.9).
Lemma 5.5. Suppose F is τ -spectrally bounded by g. Then,
kG⊤ (F F ⊤ )+ Gk ≤ 2τ .

(5.4)

Proof. Let K be the column span of F . Let Q = F F ⊤ . Multiplying (Q+ )1/2 on both sides of
equation (2.9), we obtain that
(Q+ )1/2 GG⊤ (Q+ )1/2  τ (IdK + σm (Q)Q+ )

 τ (IdK + σm (Q)kQ+ kIdK )

It follows that k(Q+ )1/2 Gk
kG⊤ (Q+ )1/2 (Q+ )1/2 Gk ≤ 2τ .

≤

√

 2τ IdK

2τ , which in turns implies that kG⊤ (F F ⊤ )+ Gk

=

We also need to bound the {1, 2}{3} norm of the following systematic error tensor. This is
important because we want to bound the spectral norm of the perturbation after the whitening
operation.
15

Lemma 5.6 (Variant of [MSS16, Theorem 6.1]). Let Γ, ∆, Θ ∈ Rd×ℓ . Let γi , δi , θi be the i-th
column of Γ, ∆, Θ, respectively. Then,
X

i∈[ℓ]

γi ⊗ δi ⊗ θi

≤ kΓk · kΘk · k∆k1→2 ≤ kΓk · kΘk · k∆k

(5.5)

{1,2}{3}

Proof of Lemma 5.6. Using the definition of k·k{1,2}{3} we have that
X
i∈[]

γi ⊗ δi ⊗ θi

=

X

(γi ⊗ δi )θi⊤

(5.6)

i∈[ℓ]

{1,2}{3}

≤

=

X

1/2

(γi ⊗ δi )(γi ⊗ δi )

X

1/2

θi θi⊤

i∈[ℓ]

i∈[ℓ]

(by Cauchy-Schwarz inequality)
1/2

X

(γi γi⊤ )
i∈[ℓ]

⊗

(δi δi⊤ )

kΘk

Next observe that we have that for any i, δi δi⊤  (max kδi k2 )Id and therefore,
(γi γi⊤ ) ⊗ (δi δi⊤ )  γi γi⊤ ⊗ (max kδi k2 )Id .

(5.7)

It follows that
X

i∈[r]

γi ⊗ δi ⊗ θi

≤
{1,2}{3}

X

i∈[r]

1/2
2

γi γi⊤ ⊗ (max kδi k )Id

kΘk

= kΓk · kΘk · k∆k1→2 .

With this in mind, we prove the main theorem:
1/2 A, B̃ = (Q+ )1/2 B, C̃ = (Q+ )1/2 C. Moreover, let Γ̃ =
Proof of Theorem 5.4. Let Ã = (Q+
a)
c
b
1/2 Γ and define ∆,
˜
(Q+
)
Θ̃
similarly.
Let
ã
,
b̃
,
c̃
,
γ̃
,
δ̃
,
θ̃
be
their
columns.
Then we have that T̃
1
i
i
i
i
i
a
as defined in Algorithm 3 satisfies

T̃ =

r
X
i=1

ãi ⊗ b̃i ⊗ c̃i +

ℓ
X
i=1

γ̃i ⊗ δ̃i ⊗ θ̃i + Ẽ

(5.8)

1/2 ⊗ (Q+ )1/2 ⊗ (Q+ )1/2 · E. We will show that T̃ meets the condition of Thewhere Ẽ = (Q+
c
a)
b
1/2 A is
orem 2.6. Since Qa is an ε-approximate whitening matrix of A, by Definition, Ã = (Q+
a)
ε-approximately orthonormal . Similarly, B̃, C̃ are ε-approximately orthonormal . √
Γ is τ -spectrally
bounded by Qa , hence by Lemma 5.5, we have that kΓ̃k ≤ 2τ . Similarly,
√
kΘk, k∆k ≤ 2τ . Applying Lemma 5.6, we have,
ℓ
X
i=1

γ̃i ⊗ δ̃i ⊗ θ̃i

{1,2}{3}

16

≤ (2τ )3/2

(5.9)

−3/2 kEk
1/2 k · k(Q+ )1/2 k · k(Q+ )1/2 kkEk
Moreover, we have kẼk{1,2}{3} ≤ k(Q+
{1,2}{3} ≤ σ
{1,2}{3} ,
c
c
a)
where σ = min{σmin (Qa ), σmin (Qb ), σmin (Qc )}.PTherefore, using Theorem 2.6 (with ai , bi , ci there
replaced by ãi , b̃i , c̃i , and Z there replaced by ℓi=1 γ̃i ⊗ δ̃i ⊗ θ̃i + Ẽ), we have that a set of vectors
{ăi , b̆i , c̆i } that are ε-close to {ãi , b̃i , c̃i } with ε = (2τ )3/2 + σ −3/2 kẼk{1,2}{3} . Therefore, we obtain
that kai − Q1/2 ăi k ≤ kQa k1/2 ε. Similarly we can control the error for bi and ci and complete the
proof.

6

Conclusions

We have presented theoretical progress on the longstanding open problem of presenting a
polynomial-time algorithm for learning noisy-or networks given sample outputs from the network.
In particular it is enouraging that linear algebraic methods like tensor decomposition can play a
role. Earlier there were no good approaches for this problem; even heuristics fail for realistic sizes
like n = 1000.
Can sample complexity be reduced, say to subcubic? (Cubic implies more than one billion
examples for networks with 1000 outputs.) Possibly this requires exploiting some hierarchichal
structure –e.g. groupings of diseases and symptoms— in practical noisy-OR networks but exploring
such possibilities using the current version of QMR-DT is difficult because it has been scrubbed of
labels for diseases and symptoms.)
Various more practical versions of our algorithm are also easy to conceive and will be tested
in the near future. This could be somewhat analogous to topic modeling, for which discovery of
provable polynomial-time algorithms soon led to very efficient algorithms.

Acknowledgments:
We thank Randolph Miller and Vanderbilt University for providing the current version of this
network for our research.

References
[AGH+ 13] Sanjeev Arora, Rong Ge, Yonatan Halpern, David M Mimno, Ankur Moitra, David
Sontag, Yichen Wu, and Michael Zhu. A practical algorithm for topic modeling with
provable guarantees. In ICML (2), pages 280–288, 2013.
[AGH+ 14] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine
Learning Research, 15(1):2773–2832, 2014.
[AGM12]

Sanjeev Arora, Rong Ge, and Ankur Moitra. Learning topic models–going beyond svd.
In Foundations of Computer Science (FOCS), 2012 IEEE 53rd Annual Symposium on,
pages 1–10. IEEE, 2012.

[AGMM15] Sanjeev Arora, Rong Ge, Tengyu Ma, and Ankur Moitra. Simple, efficient, and neural
algorithms for sparse coding. In Proceedings of The 28th Conference on Learning
Theory, pages 113–149, 2015.
[Chi96]

David Maxwell Chickering. Learning bayesian networks is np-complete. In Learning
from data, pages 121–130. Springer, 1996.

17

[DK70]

Chandler Davis and William Morton Kahan. The rotation of eigenvectors by a perturbation. iii. SIAM Journal on Numerical Analysis, 7(1):1–46, 1970.

[GHK15]

Rong Ge, Qingqing Huang, and Sham M Kakade. Learning mixtures of gaussians in
high dimensions. In Proceedings of the Forty-Seventh Annual ACM on Symposium on
Theory of Computing, pages 761–770. ACM, 2015.

[HK13]

Daniel Hsu and Sham M Kakade. Learning mixtures of spherical gaussians: moment
methods and spectral decompositions. In Proceedings of the 4th conference on Innovations in Theoretical Computer Science, pages 11–20. ACM, 2013.

[HS06]

Geoffrey E Hinton and Ruslan R Salakhutdinov. Reducing the dimensionality of data
with neural networks. Science, 313(5786):504–507, 2006.

[HS13]

Yoni Halpern and David Sontag. Unsupervised learning of noisy-or bayesian networks.
In Uncertainty in Artificial Intelligence, page 272. Citeseer, 2013.

[Ips98]

Ilse CF Ipsen. Relative perturbation results for matrix eigenvalues and singular values.
Acta numerica, 7:151–201, 1998.

[JGJS99]

Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul. An
introduction to variational methods for graphical models. Machine learning, 37(2):183–
233, 1999.

[JHS13]

Yacine Jernite, Yonatan Halpern, and David Sontag. Discovering hidden variables in
noisy-or networks using quartet tests. In Advances in Neural Information Processing
Systems, pages 2355–2363, 2013.

[Li97]

Ren-Cang Li. Relative perturbation theory. iii. more bounds on eigenvalue variation.
Linear algebra and its applications, 266:337–345, 1997.

[Li98a]

Ren-Cang Li. Relative perturbation theory: I. eigenvalue and singular value variations.
SIAM Journal on Matrix Analysis and Applications, 19(4):956–982, 1998.

[Li98b]

Ren-Cang Li. Relative perturbation theory: Ii. eigenspace and singular subspace variations. SIAM Journal on Matrix Analysis and Applications, 20(2):471–492, 1998.

[MPJM82] Randolph A Miller, Harry E Pople Jr, and Jack D Myers. Internist-i, an experimental computer-based diagnostic consultant for general internal medicine. New England
Journal of Medicine, 307(8):468–476, 1982.
[MR05]

Elchanan Mossel and Sébastien Roch. Learning nonsingular phylogenies and hidden
markov models. In Proceedings of the thirty-seventh annual ACM symposium on Theory
of computing, pages 366–375. ACM, 2005.

[MSS16]

T. Ma, J. Shi, and D. Steurer. Polynomial-time Tensor Decompositions with Sum-ofSquares. ArXiv e-prints, October 2016.

[SC91]

Michael Shwe and Gregory Cooper. An empirical analysis of likelihood-weighting simulation on a large, multiply connected medical belief network. Computers and Biomedical
Research, 24(5):453–475, 1991.

[Smo86]

Paul Smolensky. Information processing in dynamical systems: Foundations of harmony theory. Technical report, DTIC Document, 1986.
18

[Ste77]

GW Stewart. On the perturbation of pseudo-inverses, projections and linear least
squares problems. SIAM review, 19(4):634–662, 1977.

[Ste90]

Gilbert W Stewart. Matrix perturbation theory. 1990.

[Wed72]

Per-Åke Wedin. Perturbation bounds in connection with singular value decomposition.
BIT Numerical Mathematics, 12(1):99–111, 1972.

A

Formal expression for the PMI tensor

In this section we formally derive the expressions for the PMI tensors and matrices, which we only
informally did in Section 2.
As a notational convenience for l ∈ N, we will denote by P̃l the matrix which has as columns the
vectors 1 − exp(−lWk ), k ∈ [m]. Furthermore, for a subset Sa ⊆ [n], we will introduce the notation
⊤

X
X 
=
(1 − exp(−lWk )Sa ) (1 − exp(−lWk )Sa )⊤
(P̃l )k,Sa (P̃l )k,Sa
Pl,Sa =
k∈[m]

k∈[m]

These matrices will appear naturally in the expressions for the higher-order terms in the Taylor
expansion for the PMI matrix and tensor.
We first compute the formally the moments of the noisy-or model.
Lemma A.1. We have
log Pr[si = 0] =

X

k∈[m]

∀i 6= j log Pr [si = 0 ∧ sj = 0] =
∀ distinct i, j, k ∈ [n], log Pr [si = 0 ∧ sj = 0 ∧ sℓ = 0] =

X

k∈[m]

X

k∈[m]

log (1 − ρ(1 − exp(Wik )))
log (1 − ρ(1 − exp(Wik + Wjk )))
log (1 − ρ(1 − exp(Wik + Wjk + Wℓk )))

Proof of Lemma A.1. We only give the proof for the second equation. The rest can be shown
analogously.
h
i
log Pr [si = 0 ∧ sj = 0] = E [Pr[si = 0|d] · Pr[sj = 0|d]] = E exp(−(Wi + Wj )⊤ d)
Y
=
E [exp(−(Wik + Wjk )dk ))]
k∈[m]

=

Y

k∈[m]

(1 − ρ(1 − exp(−(Wik + Wjk )))) .

With this in mind, we give the expression for the PMI tensor along with all the higher-order
terms.
Proposition A.2. For any equipartition Sa , Sb , Sc of [n], the restriction of the PMI tensor
PMITSa ,Sb ,Sc satisfies, for any L ≥ 2,
PMITSa ,Sb ,Sc

L
X
X
ρ
Fk,Sa ⊗Fk,Sb ⊗Fk,Sc + (−1)l+1
=
1−ρ
k∈[m]

l=2

19

1
l



ρ
1−ρ

l ! X

k∈[m]

(P̃l )k,Sa ⊗(P̃l )k,Sb ⊗(P̃l )Sc +EL
(A.1)

where
kEL k{1,2},{3} ≤

(mn)3
L



1−

ρ
1−ρ



L

ρ
1−ρ

L

Proof. The proof will proceed by Taylor expanding the log terms. Towards that, using Lemma A.1,
we have :
PMITijl =
X
(1 − ρ(1 − exp(−Wik − Wjk ))) (1 − ρ(1 − exp(−Wik − Wlk ))) (1 − ρ(1 − exp(−Wjk − Wlk )))
log
(1 − ρ(1 − exp(−Wik − Wjk − Wlk ))) (1 − ρ(1 − exp(−Wik ))) (1 − ρ(1 − exp(−Wjk ))) (1 − ρ(1 − exp(−Wlk )))
k∈[m]

By the Taylor expansion of log(1 − x), we get that
PMIijl =
∞
X
1 X t
−
ρ (((1 − exp(−Wik − Wjk )))t + ((1 − exp(−Wik − Wlk )))t + ((1 − exp(−Wjk − Wlk )))t −
t
t=1
k∈[m]

(1 − exp(−Wik ))t − (1 − exp(−Wjk ))t − (1 − exp(−Wlk ))t − (1 − exp(−Wik − Wjk − Wlk ))t )

Furthermore, note that
((1 − exp(−Wik − Wjk )))t + ((1 − exp(−Wik − Wlk )))t + ((1 − exp(−Wjk − Wlk )))t −

(1 − exp(−Wik ))t − (1 − exp(−Wjk ))t − (1 − exp(−Wlk ))t − (1 − exp(−Wik − Wjk − Wlk ))t =
t  
X
t
(−1)l (1 − exp(−lWik )) (1 − exp(−lWjk )) (1 − exp(−lWlk ))
l
l=1

by simple regrouping of the terms. By exchanging l and t, we get

PMITSa ,Sb ,Sc

  X

t1 t
(1 − exp(−lWk ))Sa ⊗ (1 − exp(−lWk ))Sb ⊗ (1 − exp(−lWk ))Sc
ρ
=
(−1)
t l
l=1 t≥l
k∈[m]

l ! X
∞
X
1
ρ
(1 − exp(−lWk ))Sa ⊗ (1 − exp(−lWk ))Sb ⊗ (1 − exp(−lWk ))Sc
=
(−1)l+1
l 1−ρ
∞ X
X

l+1

l=1

k∈[m]

(A.2)

where the last equality holds by noting that
X 1 t  1  ρ l
ρt
=
t l
l 1−ρ
t≥l

The term corresponding to t = 1 is easily seen to be
ρ X
Fk,Sa ⊗ Fk,Sb ⊗ Fk,Sc
1−ρ
k∈[m]

therefore we to show the statement of the lemma, we only need bound the contribution of the terms
with l ≥ L.
20

Toward that, note that ∀l, kk1 − exp(−lWk )k ≤ n. Hence, we have by Lemma 5.6,
m
X
k=1

(1 − exp(−lWk ))Sa ⊗ (1 − exp(−lWk ))Sb ⊗ (1 − exp(−lWk ))Sc

{1,2},{3}

≤ (mn)3

Therefore, subadditivity of the {12}, {3} norm gives

l 
ρ
∞
X
 1−ρ  X
(1 − exp(−lWk ))Sa ⊗ (1 − exp(−lWk ))Sb ⊗ (1 − exp(−lWk ))Sc
(−1)l+1 

l
l=L

k∈[m]



∞
X

≤ (mn)

3

l=L

ρ
1−ρ

l

l 



l
∞ 
(mn)3
ρ
 (mn)3 X
=
≤
L
1−ρ
L

1−

l=L

ρ
1−ρ



L

ρ
1−ρ

{1,2},{3}

L

which gives us what we need.
A completely analogous proof gives a similar expression for the PMI matrix:
Proposition A.3. For any subsets Sa , Sb of [n], s.t. Sa ∩ Sb = ∅, the restriction of the PMI matrix
PMISa ,Sb satisfies, for any L ≥ 2,
PMISa ,Sb

L
X
X
ρ
⊤
+
(−1)l+1
Fk,Sa Fk,S
=
b
1−ρ
k∈[m]

l=2

1
l



ρ
1−ρ

where
kEL k{1,2},{3} ≤

B

(mn)2
L



1−

l ! X

(P̃l )k,Sa ((P̃l )k,Sb )⊤ + EL

(A.3)

k∈[m]

ρ
1−ρ



L

ρ
1−ρ

L

Spectral properties of the random model

The goal of this section is to prove that the random model specified in Section 1 satisfies the
incoherence property 3.2 on the weight matrix and the spectral boundedness property of the PMI
tensor. (Recall, the former is required for the whitening algorithm, and the later for the tensor
decomposition algorithm.)
Before delving into the proofs, we will need a few simple bounds on the singular values of Pl .
Lemma B.1. Let Sa ⊆ [n], s.t. |Sa | = Ω(n). With probability 1 − exp(− log2 n) over the choice of
W , and for all l = O(poly(n)),
σmin (Pl,Sa ) & np
and
σmax (Pl,Sa ) . mnp2

21

Proof. Let us proceed to the lower bound first.
If we denote by L the matrix which has as columns (P̃l )k,Sa , k ∈ [m], it’s clear that Pl,Sa = LL⊤ .
Since
σmin (LL⊤ ) = σmin (L⊤ L)
we will proceed to bound the smallest eigenvalue of L⊤ L.
Note that
X
L⊤ L =
(1 − exp(−lW k ))(1 − exp(−lW k ))⊤
k∈Sa

Since the matrices (1 − exp(−lW k ))(1 − exp(−lW k ))⊤ are independent, the bound will follow
from a matrix Bernstein bound. Denoting
i
h
Q = E (1 − exp(−lW k ))(1 − exp(−lW k ))⊤

by a simple calculation we have

h
i2

i
 h
Q = p2 E 1 − exp(−lW̃ ) 11⊤ + pE (1 − exp(−lW̃ ))2 − p2 E [1 − exp(−lw̃)]2 Idm

(B.1)

where 1 is the all-ones vector of size m, and Idm is the identity of the same size. Furthermore, W̃
is a random variable following the distribution D of all the W̃i,j .
Note that (1.2) together with the assumption ν = Ω(1) gives
= Ω(p)

Pσmin (Q)
i
i ⊤ = |S | · Id , and
i
−1/2
Let Z = Q
(1 − exp(−lWi )). Then we have that E
a
m
i∈Sa Z (Z )
with high probability, kZ i k2 ≤ 1/σmin (Q) · kF i k2 . m and it’s a sub-exponential random variable.
Moreover,
#
"
2 
X
X X
Z i (Z i )⊤ ≤ mn .
.m E
Z i (Z i )⊤
r2 =
E
i

i

Therefore, by Bernstein inequality we have that w.h.p,
X

i∈Sa

Z i (Z i )⊤ − nIdm .

It follows that

X

i∈Sa

p

r 2 log n + max kZ i k2 log n =

p

mn log n + m log n .


p

mn log n Idm  nIdm .
Z i (Z i )⊤  n − O

which in turn implies that Pl,Sa  nQ. But this immediately implies σmin (Pl,Sa ) & np with high
probability. Union bounding over all l, we get the first part of the lemma.
The upper bound will be proven by a Chernoff bound. Note that the matrices
(1 − exp(−lWk ))Sa (1 − exp(−lWk ))⊤
Sa , k ∈ [m]

2
are independent. Furthermore, k(1−exp(−lWk ))Sa (1−exp(−lWk ))⊤
Sa k ≤ pn with high probability,
2
⊤
and the variable k(1 − exp(−lWk ))Sa (1 − exp(−lWk ))Sa k is sub-exponential. Finally,
#
"m
X
2
((1 − exp(−lWk ))Sa (1 − exp(−lWk ))⊤
r2 = E
Sa )
k=1

≤

m
X
k=1

i
h
2
E ((1 − exp(−lWk )Sa )(1 − exp(−lWk ))⊤
)
Sa

h
i
2 2
≤ m k1 − exp(−lWk )Sa k2 E ((1 − exp(−lWk ))Sa (1 − exp(−lWk ))⊤
Sa ) ≤ mn p
22

Similarly as in the lower bound,

h
i2
i
 h
E[Pl,Sa ] = p2 E 1 − exp(−lW̃ ) 11⊤ + pE (1 − exp(−lW̃ ))2 − p2 E [1 − exp(−lw̃)]2 Id|Sa |

where 1 is the all-ones vector of size |Sa |, and Id|Sa | is the identity of the same size. Again, W̃ is a
random variable following the distribution D of all the W̃i,j . This immediately gives
p
Pl  E[Pl ] + r log nId|Sa |  mnp2 + mn2 p2 log nId|Sa |  O(mnp2 )Id|Sa |
A union bound over all values of l gives the statement of the Lemma.

B.1

Incoherence of matrix F

First, we proceed to show the incoherence property 3.2 on the weight matrix.
Lemma B.2. Suppose n is a multiple of 3. Let F = U ΣV be the singular value decomposition of
F . Let Sa , Sb , Sc be a uniformly random equipartition of the rows of [n]. Suppose F is µ-incoherent
with µ ≤ n/(m log n). Then, with high probability over the choice of and Sa , Sb , Sc , we have for
every i ∈ {a, b, c},
r

1
µm
Si ⊤ Si
U − Id .
log n .
U
3
n
Proof. Let S = Sa . Then, since U ⊤ U = Idm , we have
#
"
X
1
(U i )(U i )⊤ = · Idm .
E
3
i∈S

By the assumption on the row norms of U , kU i (U i )⊤ k2 = kU i k2 ≤ µ m
n . By the incoherence
assumption, we have that maxi kU i k2 ≤ µm/n.
We also note that U i ’s are negatively associated random variables. Therefore by the matrix
Chernoff inequality for negatively associated random variables, we have with high probability,
r
i
h
µm log n
S ⊤ S
S ⊤ S
.
(U ) U − E (U ) U
.
n
But, an analogous argument holds for Sb , Sc as well – so by a union bound over k, we complete
the proof.
Lemma B.3. Suppose n & m log n. Under the generative assumption in Section 1 for W , we have
that we have that F = 1 − exp(−W ) is O(1)-incoherent.

Proof. We have that F F ⊤ = U Σ2 U T and therefore, kF i k2 = Σ2i,i kU i k2 . This in turn implies that
kU i k2 ≤

1
1
Fi = 2
kF i k2 .
min Σ2ii
σmin (F )

√
2
Since F i ≤ pm+ pm ≤ 2pm with high probability, we only need to bound σmin (F ) from below.
2 (F ) = σ
⊤
⊤
Note that σmin
min (F F ). Therefore it suffices to control σmin (F F ). But by Lemma B.1
2
we have σmin (F ) & np.
Therefore, we have that
m
1
kF i k2 = O( )
kU i k2 ≤ 2
n
σmin (F )

23

B.2

Spectral boundedness

The main goal of the section is to show that the bias terms in the PMI tensor are spectrally
bounded by the PMI matrix (which we can estimate from samples). Furthermore, we show that
we can calculate an approximate whitening matrix for the leading terms of the PMI tensor using
the machinery in Section .
The main proposition we will show is the following:
Proposition
B.4.

 Let W be sampled according to the random model in Section 1 with ρ =

log n
1
. Let Sa ⊆ [n], |Sa | = Ω(n). If RSa is the matrix that has as columns
o log n , p = ω √
mn
 

l 1/3
ρ
the vectors 1l 1−ρ
(P̃l )j,Sa , l ∈ [2, L], j ∈ [m] for L = O(poly(n)), and A is the matrix that
1/3

ρ
(P̃1 )j,Sa for j ∈ [m], with high probability it holds that RSa RS⊤a
has as columns the vectors 1−ρ
is O(ρ2/3 log n)-spectrally bounded by A.

The main element for the proposition is the following Lemma:
Lemma B.5. For any set Sa ⊆ [n], |Sa | = Ω(n), with high probability over the choice of W , for
⊤ is O(log n)-spectrally bounded by P ′
every ℓ, ℓ′ = O(poly(n)) Pl,Sa Pl,S
l ,Sa .
a
Before proving the Lemma, let us see how the proposition follows from it:
Proof of B.4. We have
RSa RS⊤a

=

L
X
l=2

1
l



ρ
1−ρ

l !2/3

Pl,Sa

By Lemma B.5, we have that ∀l > 1, P̃l,Sa is τ -spectrally bounded by P̃1,Sa , for some τ =
O(log n). Hence,

RSa RS⊤a



∞
X
l=2

1
l



ρ
1−ρ

l !2/3

τ (P1,Sa + σmin (P1,Sa ))

- ρ4/3 τ (P1,Sa + σmin (P1,Sa ))
ρ 2/3
) P1,Sa , the claim of the Proposition follows.
Since AA⊤ = ( 1−ρ
It is clear analogous statements hold for Sb and Sc .

Finally, we proceed to show Lemma B.5.
For notational convenience, we will denote by Jm×n the all ones matrix with dimension m × n.
(We will omit the dimensions when clear from the context.)
This statement will immediately follow from the following two lemmas:
Lemma B.6. For any set Sa ⊆ [n], |Sa | = Ω(n), with probability 1 − exp(− log2 n) over the choice
of W , for all ℓ ≤ O(poly(n)),
5
Pl,Sa  10np log nId + mp2 J
2

24

Lemma B.7. For any set Sa ⊆ [n], |Sa | = Ω(n), with probability 1 − exp(− log2 n) over the choice
of W , ∀ℓ = poly(n),
Pl,Sa + 6np log nId - mp2 J

Before showing these lemmas, let us see how Lemma B.5 is implied by them:
Proof of Lemma B.5. Let κ be the constant in B.7, s.t. Pl,Sa + 6np log nId - mp2 J. Putting the
bounds from Lemmas B.6 and B.7 together along with a union bound, we have that with high
probability, ∀l, l′ = O(poly(n))


15
5
Pl,Sa − κPl′ ,Sa  10np log n + κnp log n Id  O(mp log n)Id
2
2
But, note that σmin (Pl′ ) = Ω(np), by Lemma B.1. Hence, Pl − 52 κPl′ ,Sa  r log nσmin (Pl′ ), for
some sufficiently large constant r. This implies
5
Pl − r log nPl′ ,Sa  Pl − κPl′  r log nσmin (Pl′ ,Sa )
2
from which the statement of the lemma follows.
We proceed to the first lemma:
Proof of Lemma B.6. To make the notation less cluttered, we will drop l and Sa and use P = Pl,Sa .
Furthermore, we will drop Sa when referring to columns of P̃ so we will denote P̃k = P̃k,Sa .
Let’s denote by e = √ 1 1 . Let’s furthermore denote Id1 = ee⊤ , and Id−1 = Id − ee⊤ . Note
|Sa |

first that trivially, since Id1 + Id−1 = Id,

P = (Id1 + Id−1 )P (Id1 + Id−1 )

(B.2)

Furthermore, it also holds that
1
1
0  (2Id−1 − Id1 )P (2Id−1 − Id1 )
2
2
1
= Id1 P Id1 + 4Id−1 P Id−1 − Id−1 P Id1 − Id1 P Id−1
4
where the first inequality holds since
1
1
(2Id−1 − Id1 )P (2Id−1 − Id1 ) =
2
2



1
(2Id−1 − Id1 )P̃
2



1
(2Id−1 − Id1 )P̃
2

⊤

From this we get that
Id−1 P Id1 + Id1 P Id−1 

1
Id1 P Id1 + 4Id−1 P Id−1
4

(B.3)

We proceed to upper bound both the terms on the RHS above. More precisely, we will show
that
Id1 P Id1 ≤ 2mp2 nJ
(B.4)
Id−1 P Id−1  2np log nId
25

(B.5)

Let us proceed to showing (B.4). The LHS can be rewritten as


Id1 P Id1 = ee⊤ e⊤ P̃ P̃ ⊤ e

Note that

m
X

1
e⊤ P̃ P̃ ⊤ e =
n

k=1

h1, P̃k i2

!

All the terms h1, P̃k i2 are independent and satisfy and E[h1, P̃k i2 ] ≤ (E[h1, P̃k i])2 ≤ p2 n2 . By
Chernoff, we have that
m
X
p
h1, P̃k i2 ≤ mp2 n2 + mp2 n2 log n ≤ 2mp2 n2
k=1

log n
).
with probability at least 1 − exp(− log2 n), where the second inequality holds because p = ω( √
mn

Hence, e⊤ P̃ P̃ ⊤ e ≤ 2mp2 n with high probability.
We proceed to (B.5), which will be shown by a Bernstein bound. Towards that, note that
E[(Id−1 P̃k (Id−1 P̃k )⊤ ] = Id−1 E[P̃k (P̃k )⊤ ]Id−1


= Id−1 p2 E[1 − exp(W̃ )]2 11⊤ + (p − p2 )E[(1 − exp(W̃ ))2 ]Id Id−1
 pId

where W̃ is a random variable following the distribution D of all the W̃i,j . The second line can be
seen to follow from the independence of the coordinates of P̃k according to our model.
Furthermore, with high probability kId−1 P̃k k22 ≤ kP̃k k22 ≤ np and the random variable
kId−1 P̃k k2 is sub-exponential Finally,
#
"m
#
"m
2
X
X
2
(Id−1 P̃k )(Id−1 Γk )P̃
≤ P̃k
r2 = E
(I−1 P̃k )(Id−1 P̃k )⊤
E
2

k=1

k=1

≤ npmkE[(Id−1 P̃k )(Id−1 P̃k )⊤ ]k ≤ np2 m

Therefore, applying a matrix Bernstein bound, we get
Id−1 P Id−1  mpId + np log nId + r log nId  2np log nId
with high probability.
Combining this with (B.2) and (B.3), we get
5
P  Id1 P Id1 + 5Id−1 P Id−1
4
5
 mp2 J + 10np log nId
2

Let us proceed to the second inequality, which essentially follows the same strategy:

26

Proof of Lemma B.7. Similarly as in the proof of Lemma B.6, for notational convenience, let’s
denote by P̃ the matrix which has column k the vector (1 − exp(lWk )).
Reusing the notation from Lemma B.7, we have that
P = (Id1 + Id−1 )P (Id1 + Id−1 )

(B.6)

and
1
1
0  ( Id1 + 2Id−1 )P ( Id1 + 2Id−1 )
2
2
1
′
= Id−1 P Id1 + Id1 P Id−1 + Id1 P Id1 + 4Id−1 P Id−1
4
for similar reasons as before. From this we get that
1
Id−1 P Id1 + Id1 P Id−1  − Id1 P Id1 − 4Id−1 P Id−1
4

(B.7)

Putting (B.6) and (B.7) together, we get that
3
P + 3Id−1 P Id−1  Id1 P Id1
4

(B.8)

We will proceed to show an upper bound Id−1 P Id−1  2np log nId on second term of the LHS.
We will do this by a Bernstein bound as before. Namely, analogously as in Lemma B.7,
Id−1 P Id−1 =

m
X
k=1

(Id−1 P̃k )(Id−1 P̃k )⊤

and E[(Id−1 P̃k )(Id−1 P̃k )⊤ ]  pId and and kId−1 P̃k k22 ≤ kP̃k k22 ≤ np are satisfied so
"m
#
2
X
(Id−1 P̃k )(Id−1 P̃k )⊤
≤ np2 m
r2 = E
k=1

Therefore, applying a matrix Bernstein bound, we get
Id−1 P Id−1  mpId + np log nId + r  2np log nId
with high probability.
Plugging this in in (B.8), we get


3
3
3
Id1 P Id1 = ee⊤ P̃ P̃ ⊤ ee⊤ = ee⊤ e⊤ P̃ P̃ ⊤ e
4
4
4
P

m
2
with the goal of applying Chernoff, we will lower
Since we have e⊤ P̃ P̃ ⊤ e = n1
k=1 hP̃k , 1i
P + 6np log nId 

bound E[h1, Ak i2 ]. More precisely, we will show E[h1, P̃k i]2 = Ω(n2 p2 ). In order to do this, we have
X
X
E[(P̃k )2j ] +
E[h1, P̃k i2 ] =
E[(P̃k )j ]E[(P̃k )j ′ ]
j6=j ′ ;j,j ′∈Sa

j∈Sa

=

X

j∈Sa

≥

X

j∈Sa

˜ )2 ] +
pE[(1 − exp(−lW
pE[(1 − exp(−W̃ )2 ] +

2 2

= Ω(n p )
27

X

j6=j ′ ;j,j ′∈Sa

X

j6=j ′ ;j,j ′∈Sa

˜ )]2
p2 E[1 − exp(−lW

p2 E[1 − exp(−W̃ )]2

where W̃ is a random variable following the distribution D of all the W̃i,j . and the last inequality
holds because of (1.2).
p
So by Chernoff, we get that e⊤ P̃ P̃ ⊤ e = n1 (Ω(mn2 p2 ) − n2 p2 m) = Ω(mnp2 ) with high probability. Altogether, this means
P + 6np log nId % mp2 J

C

Robust whitening

Algorithm 4 Obtaining whitening matrices
d
Inputs: Random partitioning Sa , Sb , Sc of [n]. Empirical PMI matrix PMI.
d×d
Outputs: Whitening matrices Qa , Qb , Qc ∈ R
1. Output

+

d Sc ,Sa ,
d Sa ,S (PMI
d S ,S )⊤ PMI
Qa = ρ−1/3 PMI
b
b c

d S ,Sc (PMI
d + )⊤ PMI
d Sa ,S ,
Qb = ρ−1/3 PMI
Sc ,Sa
b
b
⊤ d
d Sc ,Sa (PMI
d+
Qc = ρ−1/3 PMI
Sa ,Sb ) PMISb ,Sc

⊤ d
d Sa ,S (PMI
d+
In this section, we show the formula Qa = ρ−1/3 PMI
Sb ,Sc ) PMISc ,Sa computes an
b
⊤
approximation of the true whitening matrix AA , so that the error is ε-spectrally bounded by A.
We recall Theorem 2.10.

Theorem 2.10. Let n ≥ m and A, B, C ∈ Rn×m . Suppose Σab , Σbc , Σca ∈ Rn×n are of the form,
Σab = AB ⊤ + Eab , Σbc = BC ⊤ + Ebc , and Σca = CA⊤ + Eca .
where Eab , Ebc , Eca are ε-spectrally bounded by (A, B), (B, C), (C, A) respectively. Then, the matrix
matrix
+
Qa = Σab [Σ⊤
bc ]m Σca
⊤
+
is a good approximation of AA⊤ in the sense that Qa = Σab [Σ⊤
bc ]m Σca − AA is O(ε)-spectrally
bounded by A. Here [Σ]m denotes the best rank-m approximation of Σ.

Towards proving Theorem 2.10, an intermediate step is to understand the how the space of
singular vectors of BC ⊤ are aligned with the noisy version Σbc . The following explicitly represent
BC ⊤ + E as the form B ′ R(C ′ )⊤ + ∆′ . Here the crucial benefit to do so is that the resulting ∆′ is
small in every direction. In other words, we started with a relative error guarantees on E and the
Lemma below converts to it an absolute error guarantees on ∆′ (though the signal term changes
slightly).
Lemma C.1. Suppose B, C are n × m matrices with n ≥ m. Suppose a matrix E is ε-spectrally
bounded by (B, C), then BC ⊤ + E can be written as
BC ⊤ + E = (B + ∆B )RBC (C + ∆C )⊤ + ∆′BC ,

28

where ∆B , ∆C , ∆′BC are small and RBC is close to identity in the sense that,
k∆B k ≤ O(εσmin (B))
k∆C k ≤ O(εσmin (C))

k∆′BC k ≤ O(εσmin (B)σmin (C))

kRBC − Idk ≤ O(ε)

Proof. The key intuition is if the perturbation is happening in the span of columns of B and C,
they cannot change the subspace. By Definition 2.9, we can write E as
⊤
E = B∆1 C ⊤ + B∆⊤
2 + ∆3 C + ∆4 .

Now since k∆1 k ≤ ε < 1, we know (Id − ∆1 ) is invertible, so we can write
(BC ⊤ + E) = (B + ∆3 (Id + ∆1 )−1 )(Id + ∆1 )(C + ∆2 (Id − ∆1 )−⊤ )⊤ − ∆3 (Id − ∆1 )−1 ∆⊤
2 + ∆4 .
This is already in the desired form as we can let ∆B = ∆3 (Id + ∆1 )−1 , RBC = (Id + ∆1 ),
∆C = ∆2 (Id − ∆1 )−⊤ , and ∆′BC = −∆3 (Id − ∆1 )−1 ∆⊤
2 + ∆4 . By Weyl’s Theorem we know
−1
ε
σmin (B). Other terms can be
σmin (Id + ∆1 ) ≥ 1 − ε, therefore k∆B k ≤ k∆3 kσmin
(Id + ∆) ≤ 1−ε
bounded similarly.
Now we prove that the top m approximation of BC ⊤ + E has similar column/row spaces as
BC ⊤ . Let UB be the column span of B, UB′ be the column span of (B + ∆B ), and UB′′ be the top m
left singular subspace of (BC ⊤ + E). Similarly we can define UC , UC′ , UC′′ to be the column spans
of C, C + ∆C and the top m right singular subspace of (BC ⊤ + E).
For B + ∆B , we can apply Weyl’s Theorem and Wedin’s Theorem. By Weyl’s Theorem we
know σmin (B + ∆B ) ≥ σmin (B) − k∆B k ≥ (1 − O(ε))σmin (B). By Wedin’s Theorem we know UB′
is O(ε)-close to UB . Similar results apply to C + ∆C .
Now we know σmin ((B + ∆B )RBC (C + ∆C )⊤ )) ≥ σmin (B + ∆B )σmin (RBC )σmin (C + ∆C ) ≥
Ω(σmin (B)σmin (C). Therefore we can again apply Wedin’s Theorem, considering (B+∆B )RBC (C+
∆C )⊤ ) as the original matrix and ∆′BC as the perturbation. As a result, we know UB′′ is O(ε) close
to UB′ , UC′′ is O(ε) close to UB′ . The distance between UB , UB′′ (and UC , UC′′ ) then follows from
triangle inequality.
As a direct corollary of Lemma C.1, we obtain that the BC ⊤ and BC ⊤ + E have similar
subspaces of singular vectors.
Corollary C.2. In the setting of Lemma C.1, let [BC ⊤ + E]m be the best rank-m approximation
of BC ⊤ + E. Then, the span of columns of [BC ⊤ + E]m is O(ε)-close to the span of columns of
B, span of rows of [BC ⊤ + E]m is O(ε)-close to the span of columns of C.
Furthermore, we can write [BC ⊤ + E]m = (B + ∆B )RBC (C + ∆C )⊤ + ∆BC . Here ∆B , ∆C and
RBC as defined in Lemma C.1 and ∆BC satisfies k∆BC k ≤ O(εσmin (B)σmin (C)).
Proof. Since [BC ⊤ + E]m is the best rank-m approximation, because (B + ∆B )RBC (C + ∆C )⊤ ) is
a rank m matrix, in particular we have
kBC ⊤ + E − [BC ⊤ + E]m k ≤ kBC ⊤ − (B + ∆B )RBC (C + ∆C )⊤ )kk∆′BC k.

29

Therefore
k∆BC k = k[BC ⊤ + E]m − (B + ∆B )RBC (C + ∆C )⊤ )k

≤ kBC ⊤ + E − [BC ⊤ + E]m k + kBC ⊤ + E − (B + ∆B )RBC (C + ∆C )⊤ )k

≤ 2k∆′BC k.

+
In order to fix this problem, we notice that the matrix [Σ⊤
bc ]m is multiplied by Σab on the left
⊤
⊤
+
and Σca on the right. Assuming Σab = AB , Σca = CA , we should expect [Σ⊤
bc ]m to “cancel” with
the B ⊤ factor on the left and the C factor on the right, giving us AA⊤ . Therefore, we should really
⊤
+
measure the error of the middle term [Σ⊤
bc ]m after left multiplying with B and right multiplying
with C. We formalize this in the following lemma:
+
⊤ +
Lemma C.3. Suppose Σbc is as defined in Theorem 2.10, let ∆ = [Σ⊤
bc ]m − [CB ] , then we have

kB ⊤ ∆Ck = O(ε),
k∆Ck ≤ O(

ε
),
σmin (B)

ε
),
σmin (C)
ε
k∆k ≤ O(
).
σmin (B)σmin (C)

kB ⊤ ∆k ≤ O(

We will first prove Theorem 2.10 assuming Lemma C.3.
Proof of Theorem 2.10. By Lemma C.1, we know Σab can be written as
(A + ∆1A )RAB (B + ∆1B )⊤ + ∆AB .
Similarly Σca can be written as
(C + ∆3C )RCA (A + ∆3A )⊤ + ∆CA .
Here the ∆ terms and R terms are bounded as in Lemma C.1.
+
Now let us write the matrix Σab [Σ⊤
bc ]m Σca as




(A + ∆1A )RAB (B + ∆1B )⊤ + ∆AB ([CB ⊤ ]+ + ∆BC ) (C + ∆3C )RCA (A + ∆3A )⊤ + ∆CA

+
We can now view Σab [Σ⊤
bc ]m Σca as the product of three terms, each term is the sum of two
matrices. Therefore we can expand the product into 8 terms. In each of the three pairs, we will
call the first matrix the main matrix, and the second matrix the perturbation.
In the remaining proof, we will do calculations to show the product of the main terms is close
to AA⊤ , and all the other 7 terms are small.
Before doing that, we first prove several Claims about PSD matrices

Claim C.4. If k∆k ≤ ε, then A∆A⊤  εAA⊤ . If kΓk ≤ εσmin (A), then
2 (A)Id .
εAA⊤ + εσmin

1
⊤
2 (AΓ

+ ΓA⊤ ) 

Proof. Both inequalities can be proved by consider the quadratic form. We know for any x,
x⊤ A∆A⊤ x ≤ k∆kkA⊤ xk2 ≤ εx⊤ AA⊤ x, so the first part is true.
For the second part, for any x we can apply Cauchy-Schwartz inequality
√
1
2
x⊤ (AΓ⊤ + ΓA⊤ )x = h εA⊤ x, ε−1/2 Γ⊤ xi ≤ εkA⊤ xk2 + ε−1 kΓ⊤ xk2 = x⊤ (εAA⊤ + εσmin
(A)Id)x.
2

30

Now, we will first prove the product of three main matrices is close to AA⊤ :


Claim C.5. We have (A + ∆1A )RAB (B + ∆1B )⊤ (CB ⊤ )+ (C + ∆3C )RCA (A + ∆3A )⊤ = AA⊤ +
EA , where EA is O(ε)-spectrally bounded by AA⊤ .
Proof. We will first prove the middle part of the matrix (B +∆1B )⊤ (CB ⊤ )+ (C +∆3C ) is O(ε) close to
identity matrix Id. Here we observe that both B, C have full column rank so (CB ⊤ )+ = (B ⊤ )+ C + .
Therefore we can rewrite the product as (Id + B + ∆1B )⊤ (Id + C + ∆3C ). Since k∆1B k ≤ O(εσmin (B))
by Lemma C.1 (and similarly for C), we know kB + ∆1B k ≤ O(ε). Therefore the middle part is O(ε)
bAB = RAB (B + ∆1 )⊤ (CB ⊤ )+ (C + ∆3 )RCA is O(ε)-close
close to Id. Now since ε ≪ 1 we know R
B
C
to Id.
bAB (A + ∆3 )⊤ , for this matrix we know
Now we are left with (A + ∆1A )R
A
bAB (A + ∆3A )⊤ − AA⊤ = A(R
bAB − Id)A⊤ + ∆1A R
bAB A⊤ + AR
bAB (∆3A )⊤ + ∆1A R
bAB (∆3A )⊤ .
(A + ∆1A )R

3 ⊤ 
bAB − Id)A⊤  O(ε)AA⊤ (Claim C.4); the fourth term ∆1 R
b
The first term A(R
A AB (∆A )
2 (A))Id (by the norm bounds of ∆1 and ∆3 . For the cross terms, we can bound them
O(εσmin
A
A
using the second part of Claim C.4.

Next we will try to prove the remaining 7 terms are small. We partition them into three types
depending on how many ∆ factors they have. We proceed to bound them in each of these cases.
For the terms with only one ∆, we claim:
Claim C.6. The three terms ∆AB (CB ⊤ )+ (C + ∆3C )RCA (A + ∆3A )⊤ , (A + ∆1A )RAB (B +
∆1B )⊤ ∆AB (C + ∆3C )RCA (A + ∆3A )⊤ , (A + ∆1A )RAB (B + ∆1B )⊤ (CB ⊤ )+ ∆CA are all O(ε) spectrally bounded by AA⊤ .
Proof. For the first term, note that both B, C have full column rank, and hence (CB ⊤ )+ =
(B ⊤ )+ C + . Therefore the first term can be rewritten as
[∆AB (B ⊤ )+ ][(Id + C + ∆3C )RCA ](A + ∆3A )⊤ .
By Lemma C.1, we have spectral norm bounds for ∆AB , ∆3C , ∆3A , RCA . Therefore we know
k∆AB (B ⊤ )+ k ≤ O(εσmin (A)) and [(Id + C + ∆3C )RCA ] is O(ε) close to Id.
Therefore
2 (A)) is trivially O(ε) spectrally bounded, and
k[∆AB (B ⊤ )+ ][(Id + C + ∆3C )RCA ](∆3A )⊤ k ≤ O(εσmin
[∆AB (B ⊤ )+ ][(Id + C + ∆3C )RCA ]A⊤ is O(ε) spectrally bounded by Claim C.4. The third term is
exactly symmetric.
b BC = (B+∆1 )⊤ ∆BC (C +
For the second part, we will first prove the middle part of the matrix ∆
B
∆3C ) has spectral norm O(ε). This can be done y expanding it to the sum of 4 terms, and use
appropriate spectral norm bounds on ∆BC and its products with B ⊤ and C from Lemma C.3.
b BC RCA (A + ∆3 )⊤ is O(ε) spectrally bounded by the first part
Now we can show (A + ∆1A )RAB ∆
A
of Claim C.4.
Next we try to bound the terms with two ∆ factors.

Claim C.7. The three terms
∆AB ∆BC (C + ∆3C )RCA (A + ∆3A )⊤ , ∆AB (CB ⊤ )+ ∆CA ,

(A + ∆1A )RAB (B + ∆1B )⊤ ∆BC ∆CA are all O(ε2 ) spectrally bounded by AA⊤ .

Proof. For the first term, notice that k∆BC (C + ∆3C )k is bounded by O(ε/σmin (B)) by
Lemma C.3, and k∆AB k = O(εσmin (A)σmin (B)). Therefore we know k∆AB ∆BC (C + ∆3C )RCA k ≤
O(ε2 σmin (A)), so by Claim C.4 we know this term is O(ε2 ) spectrally bounded by AA⊤ . Third
term is symmetric.
31

2 (A)),
For the second term, by Lemma C.1 we can directly bound its spectral norm by O(ε2 σmin
so it is trivially O(ε2 ) spectrally bounded by AA⊤ .

Finally, for the product ∆AB ∆BC ∆CA , we can get the spectral norm for the three factors by
2 (A)) which is trivially O(ε3 )
Lemma C.1 and Lemma C.3. As a result k∆AB ∆BC ∆CA k ≤ O(ε3 σmin
⊤
spectrally bounded by AA .
Combining the bound for all of the terms we get the theorem.
With that, we now try to prove Lemma C.3
We first prove a simpler version where the perturbation is simply bounded in spectral norm
Lemma C.8. Suppose B, C are n × m matrices and n ≥ m. Let R be an n × n matrix such that
kR − Idk ≤ ε, and E is a perturbation matrix with kEk ≤ εσmin (B)σmin (C) and (CRB ⊤ + E) is
also of rank m.
Now let ∆ = (CRB ⊤ + E)+ − (CRB ⊤ )+ , then when ε ≪ 1 we have
kB ⊤ ∆Ck = O(ε),
k∆Ck ≤ O(

ε
),
σmin (B)

ε
),
σmin (C)
ε
k∆k ≤ O(
).
σmin (B)σmin (C)

kB ⊤ ∆k ≤ O(

Proof. We first give the proof for kB ⊤ ∆Ck. Other terms are similar.
Let UB be the column span of B, and UB′ be the row span of (CRB ⊤ + E). Similarly let UC
be the column span of C and UC′ be the column span of (CRB ⊤ + E). By Wedin’s theorem, we
know UB′ is O(ε) close to UB and UC′ is O(ε) close to UC . As a result, suppose the SVD of B is
UB DB VB⊤ , we know
σmin (B ⊤ UB′ ) = σmin (VB DB UB⊤ UB′ ) ≥ (1 − O(ε))σmin (B).
The same is true for C: σmin (C ⊤ UC′ ) ≥ (1 − O(ε))σmin (C).
By the property of pseudoinverse, the column span of (CRB ⊤ + E)+ is UB′ , and the row span
of (CRB ⊤ + E)+ is UC′ , further, (CRB ⊤ + E)+ = UB′ [(UC′ )⊤ (CRB ⊤ + E)UB′ ]−1 UC′ , therefore we
can write
B ⊤ (CRB ⊤ + E)+ C = B ⊤ UB′ [(UC′ )⊤ (CRB ⊤ + E)UB′ ]−1 (UC′ )⊤ C.
Note that now the three matrices are all n×n and invertible! We can write B ⊤ UB′ = ((B ⊤ UB′ )−1 )−1
(and do the same thing for (UC′ )⊤ C. Using the fact that P −1 Q−1 = (QP )−1 , we have
B ⊤ (CRB ⊤ + E)+ C = (((UC′ )⊤ C)−1 (UC′ )⊤ (CRB ⊤ + E)UB′ (B ⊤ UB′ )−1 )−1
= (R + ((UC′ )⊤ C)−1 (UC′ )⊤ EUB′ (B ⊤ UB′ )−1 )−1 =: (R + X)−1 .
Here we defined X = ((UC′ )⊤ C)−1 (UC′ )⊤ EUB′ (B ⊤ UB′ )−1 . The spectral norm of X can be
bounded by
kXk ≤ k((UC′ )⊤ C)−1 kkEkk(B ⊤ UB′ )−1 k
−1
−1
′
= kEkσmin
(B ⊤ UB′ )σmin
(C ⊤ CB
)

≤ O(ε).

We can write B ⊤ ∆C = B ⊤ (CRB ⊤ + E)+ C − Id = (Id + (R − Id + X))−1 − Id, and we now know
k(R − Id + X)k ≤ O(ε), as a result kB ⊤ ∆Ck ≤ O(ε) as desired.
32

For the term kB ⊤ ∆k, by the same argument we have
B ⊤ (CRB ⊤ + E)+ = ((UC′ )⊤ (CRB ⊤ + E)UB′ (B ⊤ UB′ )−1 )−1 (UC′ )⊤
= ((UC′ )⊤ CR + (UC′ )⊤ EUB′ (B ⊤ UB′ )−1 )−1 (UC′ )⊤
= ((UC′ )⊤ C(R + X))−1 (UC′ )⊤
= (R + X)−1 ((UC′ )⊤ C)−1 (UC′ )⊤ .
On the other hand, we know B ⊤ (CRB ⊤ )+ = R−1 C + = R−1 ((UC )⊤ C)−1 UC⊤ . We can match
the three factors:
kR−1 − (R + X)−1 k ≤ O(ε),

k((UC )⊤ C)−1 − ((UC′ )⊤ C)−1 k ≤ O(ε/σmin (C)),
kUC − UC′ k ≤ O(ε),

kR−1 k ≤ 1 + O(ε)

k((UC )⊤ C)−1 k = O(1/σmin (C))

kUC k = 1.

Here, first and third bound are proven before. The second bound comes if we consider the SVD of
C = UC DC VC⊤ and notice that k(UC′ )⊤ UC − Idk ≤ O(ε). We can write ∆1 = R−1 − (R + X)−1 ,
∆2 = ((UC )⊤ C)−1 − ((UC′ )⊤ C)−1 , ∆3 = UC − UC′ , then we have
B ⊤ ∆ = B ⊤ (CRB ⊤ + E)+ − B ⊤ (CRB ⊤ )+

= (R−1 − ∆1 )(((UC )⊤ C)−1 − ∆2 )(UC − ∆3 )⊤ − R−1 ((UC )⊤ C)−1 UC⊤ .

Expanding the last equation, we get 7 terms and all of them can be bounded by O(ε/σmin (C)).
The bounds on k∆Ck and k∆k can be proved using similar techniques.
Finally we are ready to prove the main Lemma C.3:
⊤ , we can write the matrix before pseudoinverse
Proof of Lemma C.3. Using Lemma C.1, let E = Ebc
as
[CB ⊤ + E]m = (C + ∆C )RBC (B + ∆B )⊤ + ∆BC .

We can then apply Lemma C.8 on (C + ∆C )RBC (B + ∆B )⊤ + ∆BC . As a result, we know if we
⊤ +
let ∆′ = [CB ⊤ + E]+
m − ((C + ∆C )RBC (B + ∆B ) ) , we have the desired bound if we left multiply
⊤
with (B + ∆B ) or right multiply with (C + ∆C ).
We will now show how to prove the first bound, all the other bounds can be proved using the
same strategy:
First, we can write
′
⊤ ′
B ⊤ ∆′ C = −(B + ∆B )⊤ ∆′ (C + ∆C ) + (B + ∆B )⊤ ∆′ C + ∆⊤
B ∆ (C + ∆C ) − ∆B ∆ ∆C .

All the four terms on the RHS can be bounded by Lemma C.8 so we know kB ⊤ ∆′ Ck ≤ O(ε).
On the other hand, let ∆′′ = ((C + ∆C )RBC (B + ∆B )⊤ )+ − (C ⊤ B)+ = ∆ − ∆′ . We will prove
kB ⊤ ∆′′ Ck ≤ O(ε) and then the bound on kB ⊤ ∆Ck follows from triangle inequality.
For B ⊤ ∆′′ C, we know it is equal to
−1
(C + ∆C )+ C − Id
B ⊤ [(B + ∆B )⊤ ]+ RAB
−1
(C + ∆C )+ C − Idk ≤ O(ε)
Claim C.9. kB ⊤ [(B + ∆B )⊤ ]+ RAB

33

−1
this follows
Proof. We will show all three factors in the first term are O(ε) close to Id. For RAB
+
immediately from Lemma C.1. For (C + ∆C ) C, we know

(C + ∆C )+ C − Id = −(C + ∆C )+ ∆C .
−1
Therefore its spectral norm bound is bounded by k∆C kσmin
(C + ∆C ) = O(ε) (where the bound
on k∆C k comes from Lemma C.1).

With the claim we have now proven kB ⊤ ∆′′ Ck ≤ O(ε), therefore
kB ⊤ ∆Ck ≤ kB ⊤ (∆′ + ∆′′ )Ck ≤ kB ⊤ ∆′ Ck + kB ⊤ ∆′′ Ck ≤ O(ε).

C.1

Spectrally Boundedness and Incoherence

Here we will show under mild incoherence conditions (defined below), if an error matrix E is
ε-spectrally bounded by F F ⊤ , then the partial matrices satisfy the requirement of Theorem 2.10.
q
Theorem C.10. If F is µ-incoherent for µ ≤ n/m log2 n, then when n ≥ Ω(m log2 m), with high
probability over the random partition of F into A, B, C, we know σmin (A) ≥ σmin (F )/3 (same is
true for B, C).
As a corollary, if E is ε-spectrally bounded by F . Let a, b, c be the subsets corresponding to
A, B, C, and let Ea,b be the submatrix of E whose rows are in set a and columns are in set b. Then
Ea,b (also Eb,c , Ec,a ) is O(ε)-spectrally bounded by the corresponding asymmetric matrices AB ⊤
(BC ⊤ , CA⊤ ).
Proof. Consider the singular value decomposition of F : F = U DV ⊤ . Here U is a n × m matrix
whose columns are orthonormal, V is an m × m orthonormal matrix and D is a diagonal matrix
whose smallest diagonal entry is σmin (F ).
Consider the following way of partitioning the matrix: for each row of F , we put it into A, B
or C with probability 1/3 independently.
Now, let Xi = 1 if row i is in the matrix A, and 0 otherwise. Then Xi ’s are Bernoulli random
variables with probability 1/2. Suppose S is the set of rows in A, let UA be U restricted to rows in
A, then we have A = UA DV ⊤ . We will show with high probability
σmin (A) ≥ 1/3.
P
The key observation here is the expectation of UA⊤ UA = ni=1 Xi Ui Ui⊤ , where Ui is the i-th row
of U (represented as a column vector). Since Xi ’s are Bernoulli random variables, we know
n
n
X
1
1X
Ui Ui⊤ = Id.
Xi Ui Ui⊤ ] =
E[UA⊤ UA ] = E[
3
3
i=1

i=1

Therefore we can hope to use matrix concentration to prove that UA⊤ UA is close to its expectation.
Let Mi = Xi Ui Ui⊤ − 1/3Ui Ui⊤ . Clearly E[Mi ] = 0. By the Incoherence assumption, we know
kUi k ≤ 1/ log n. Therefore we know kMi k ≤ O(1/ log n). Also, we can bound the variance
n
n
n
X
X
X
⊤
⊤
2
⊤
Ui Ui⊤ k ≤ O(1/ log 2 n).
Xi Ui Ui Ui Ui ]k ≤ max kUi k k
Mi Mi ]k ≤ kE[
kE[
i=1

i=1

i=1

Here the last inequality is because

Pn

⊤
⊤
i=1 Xi Ui Ui Ui Ui

34

 kUi k2 Ui Ui⊤ .

Therefore by Matrix Bernstein’s inequality we know with high probability k
When this happens we know
kUA⊤ UA k

≥

σmin (E[UA⊤ UA ])

−k

n
X
i=1

Pn

i=1 Mi k

≤ 1/6.

Mi k ≥ 1/6.

p
Hence we have σmin (UA ) ≥ 1/6 > 1/3, and σmin (A) ≥ σmin (UA )σmin (D) ≥ σmin (F )/3. Note
that matrices B, C have exactly the same distribution as A so the bounds for B,C follows from
union bound.
For the corollary, if a matrix E is ε spectrally bounded, we can write it as F ∆1 F ⊤ + F ∆⊤
2 +
2 (F ). This can be done by
∆2 F ⊤ + ∆4 , where k∆1 k ≤ ε, k∆2 k ≤ εσmin (F ) and k∆4 k ≤ εσmin
considering different projections of E: let U be the span of columns of F , then F ∆1 F ⊤ term
⊤
corresponds to ProjU E ProjU ; F ∆⊤
2 term corresponds to ProjU E ProjU ⊥ ; ∆2 F term corresponds
to ProjU ⊥ E ProjU ; ∆4 term corresponds to ProjU ⊥ E ProjU ⊥ . The spectral bounds are necessary
for E to be spectrally bounded.
⊤
Now for Ea,b , we can write it as A∆1 B ⊤ + A(∆2 )⊤
b + (∆2 )a B + (∆4 )a,b , where we also take
the corresponding submatrices of ∆’s. Since the spectral norm of a submatrix can only be smaller,
we know k∆1 k ≤ ε, k(∆2 )b k ≤ εσmin (F ) ≤ 3εσmin (B), k(∆2 )a k ≤ εσmin (F ) ≤ 3εσmin (A) and
2 (F ) ≤ 9εσ
k(∆2 )a,b k ≤ εσmin
min (A)σmin (B). Therefore by Definition 2.9 we know Ea,b is 9ε
⊤
spectrally bounded by AB .

D

Proof of Theorem 3.1 and Theorem 3.3

In this section, we provide the full proof of Theorem 3.1. We start with a simple technical Lemma.
Lemma D.1. If Q is an ε-approximate whitening matrix for A, then kQk ≤
1
⊤
1−ε kAA k

1
⊤
1−ε kAA k,

σmin (Q) ≥

Proof. By the definition of approximate-whitening, we have
1 − ε ≤ σmin ((Q+ )1/2 A⊤ A(Q+ )1/2 ), σmax ((Q+ )1/2 A⊤ A(Q+ )1/2 ) ≤ 1 + ε
which implies that
1 − ε ≤ σmin ((Q+ )1/2 AA⊤ (Q+ )1/2 ), σmax ((Q+ )1/2 AA⊤ (Q+ )1/2 ) ≤ 1 + ε
by virtue of the fact that (Q+ )1/2 AA⊤ (Q+ )1/2 =
semidefinite-order notation, we get that

(Q+ )1/2 A⊤ A(Q+ )1/2

⊤

.

Rewriting in

(1 − ε)Id  (Q+ )1/2 AA⊤ (Q+ )1/2  (1 + ε)Id
Multiplying on the left and right by Q1/2 , we get
(1 − ε)Q  AA⊤  (1 + ε)Q
This directly implies

1
⊤
1+ε AA

Q

1
⊤
1−ε AA

which is equivalent to the statement of the lemma.

Towards proving Theorem 3.1, we will first prove the following proposition, which shows that
we recover the exp(−W ) matrix correctly:
35

Proposition D.2 (Recovery of exp(−W )). Under the random generative model defined in Section 1, if the number of samples satisfies
N = poly(n, 1/p/, 1/ρ)
√
the vectors W̃i , i ∈ [m] in Algorithm 1 are O(η np)-close to exp(−Wi ) where
√
η = Õ ( mpρ)

Proof. The proof will consist of checking the conditions for Algorithms 4 and 3 to work, so that we
can apply Theorems 5.4 and 2.10.
To get a handle on the PMI tensor, by Proposition A.2, for any equipartition Sa , Sb , Sc of [n],
we can it as
PMITSa ,Sb ,Sc =
L
X
ρ X
(−1)l+1
Fk,Sa ⊗ Fk,Sb ⊗ Fk,Sc +
1−ρ
k∈[m]

l=2

1
l



ρ
1−ρ

l ! X

(P̃l )k,Sa ⊗ (P̃l )k,Sb ⊗ (P̃l )k,Sc + EL

k∈[m]

1 1
We can choose L = poly(log(n, , )) to ensure
ρ p


√
kEL k = o p5/2 ρ7/3 mn2

(D.1)

Having an explicit form for the tensor, we proceed to check the spectral boundedness condition
for Algorithm 4.
Let Sa , Sb , Sc be a random equipartition. Let RSa be the matrix that has as columns the vectors

l !1/3
ρ
1
(P̃l )j,Sa , for all l ∈ [2, L], j ∈ [m] and let A be the matrix that has as columns the
l 1−ρ
1/3

ρ
Fj,Sa , for all j ∈ [m]. Since L is polynomially bounded in n, by Proposition B.4
vectors 1−ρ

we have that with high probability RSa is τ -spectrally bounded by A, for a τ = O(ρ2/3 log n).
Analogous statements hold for Sb , Sc .
Next, we verify the conditions for calculating approximate whitening matrices (Algorithm 4).
Towards applying Theorem C.10, note that if R[n] is the matrix that has as columns the vectors

l !1/3
ρ
1
(P̃l )j , for all l ∈ [2, L], j ∈ [m], and D is the matrix that has as columns the
l 1−ρ
1/3

ρ
Fj , for all j ∈ [m], R[n] then is τ -spectrally bounded by D for τ = O(ρ2/3 log n).
vectors 1−ρ
Furthermore, the matrix F is O(1)-incoherent with high probability by Lemma B.1. Hence, we can
apply Theorem C.10, the output of Algorithm 4 are matrices Qa , Qb , Qc which are τ -approximate
whitening matrices for A, B, C respectively.
Next, we will need bounds on
min(σmin (Qa ), σmin (Qb ), σmin (Qc )), max(σmax (Qa ), σmax (Qb ), σmax (Qc ))
to plug in the guarantee of Algorithm 4.

36

By Lemma D.1, we have
σmax (Qa ) ≤

1
1
kAA⊤ k . (1 + τ )kAA⊤ k, σmin (Qa ) ≥
σmin (AA⊤ ) & (1 − τ )σmin (AA⊤ )
1−τ
1+τ

However, for the random model, applying Lemma B.1,
⊤

σmin (AA ) ≥



ρ
1−ρ

2/3

np & ρ

2/3

⊤

np, σmax (AA ) ≤



ρ
1−ρ

2/3

mnp2 . ρ2/3 mnp2

(D.2)

Analogous statements hold for B and C.
Finally, we bound the error due to empirical estimates. Since ρpm = o(1),
Pr[si = 0 ∧ sj = 0 ∧ sk = 0] ≥ 1 − Pr[si = 1] − Pr[sj = 1] − Pr[sk = 1] ≥ 1 − 3pmρ = Ω(1)
Hence, by Corollary E.3, with a number of samples as stated in the theorem,
√
ˆ Sa ,S ,Sc − PMITSa ,S ,Sc k{1,2},{3} . p5/2 ρ5/2 mn2
kPMIT
b
b

(D.3)

as well.
With that, invoking Theorem 5.4 (with kEk{1,2},{3} taking into account both the EL term above,
and the above error due to sampling), the output of Algorithm 3 will produce vectors vi , i ∈ [m],
1/3

ρ
(1 − exp(−Wi )), for
s.t. vi is O(η ′ )-close to 1−ρ



ˆ Sa ,S ,Sc − PMITSa ,S ,Sc k{1,2},{3} )
η ′ . max(kQa k , kQb k , kQc k)1/2 · τ 3/2 + σ −3/2 (kEL k{1,2},{3} + kPMIT
b
b
where σ = min(σmin (Qa ), σmin (Qb ), σmin (Qc )).
Plugging in the estimates from (D.1), (D.2), (D.3) as well as τ = O(ρ2/3 log n), we get:
q
√
1/2 3/2
max(kQa k , kQb k , kQc k) τ
. mnp2 ρ2/3 (ρ2/3 log n)3/2 = ρ4/3 mnp log3/2 n
√
1 3/2
) kEL k{1,2},{3} . ρ4/3 mnp log3/2 n
σ −3/2 kEL k{1,2},{3} . (
ρnp
√
−3/2
ˆ Sa ,S ,Sc − PMITSa ,S ,Sc k{1,2},{3} . ρ4/3 mnp log3/2 n
σ
kPMIT
b

b

1/3

ρ
(1 − exp(Wi )), for all i ∈ [m].
which implies the vectors (âi , b̂i , ĉi ) are O(η ′ )-close to 1−ρ

1/3
However, this directly implies that 1−ρ
(âi , b̂i , ĉi ) are O(η ′ /ρ1/3 )-close to 1 − exp(Wi ),
ρ

i ∈ [m], which in turn implies (ãi , b̃i , c̃i ) are O(η ′ /ρ1/3 ) close to exp(Wi ).
This implies the statement of the Lemma.

Given that, we prove the main theorem. The main issue will be to ensure that taking log of the
values
√
Proof of Theorem 3.1. By Proposition D.2, the vectors Yi , i ∈ [m] in Algorithm 1 are O(η np)
√
close to exp(−Wi ) with η = Õ mpρ . Let (Yi′ )j = (Yi )j if (Yi )j > exp(−νu ) and otherwise
(Yi )′j = exp(−νu ).
Then we have that kYi′ − Wi k ≤ kYi − Wi k. By the Lipschitzness of log(·) in the region [νi , ∞]
37

we have that

It follows that

ci )j − (Wi )j | = | log(Y ′ )j − log(Wi )j | . |(Yi )′ − (Wi )j |
|(W
i
j

ci − Wi k = klog Yi′ − log Wi k . kYi′ − Wi k
kW
√
Therefore recalling kYi′ − Wi k ≤ kYi − Wi k ≤ O(η np) we complete the proof.
Proof of Theorem 3.3. The proof will follow the same outline as the proof of Theorem 3.1. The
difference is that since we only have a guarantee on the spectral boundedness of the second and
third-order term, we will need to bound the higher-order terms in a different manner. Given that
we have no information on them in this scenario, we will simply bound them in the obvious manner.
We proceed to formalize this.
The sample complexity is polynomial for the same reasons as in the proof of Theorem 3.1, so
we will not worry about it here.
We only need to check the conditions for Algorithms 4 and 3 to work, so that we can apply
Theorems 5.4 and 2.10.
Towards that, first we claim that we can write the PMI tensor for any equipartition Sa , Sb , Sc
of [n] as
PMITSa ,Sb ,Sc =


2 ! X
1
ρ
ρ X
Fk,Sa ⊗ Fk,Sb ⊗ Fk,Sc −
Gk,Sa ⊗ Gk,Sb ⊗ Gk,Sc
1−ρ
2 1−ρ
k∈[m]
k∈[m]

3 ! X
ρ
1
Hk,Sa ⊗ Hk,Sb ⊗ Hk,Sc + E
+
3 1−ρ

(D.4)

k∈[m]

where kEk{1,2},{3} ≤ ρ4 m(np)3/2 . Towards achieving this, first we claim that Proposition 5.6
implies that for any subsets Sa , Sb , Sc ,
m
X
k=1

(1 − exp(−lWk ))Sa ⊗ (1 − exp(−lWk ))Sb ⊗ (1 − exp(−lWk ))Sc

{1,2},{3}

≤ m(np)3/2

Indeed, if we put γk = (1 − exp(−lWk ))Sa , δk = (1 − exp(−lWk ))Sb , θk = (1 − exp(−lWk ))Sc ,
P
√
then we have k k γk γk⊤ k ≤ mnp, and similarly for δk . Since maxk kθk k ≤ (np)1/2 , the claim
immediately follows. Hence, (D.4) follows.
 
l 1/3
ρ
1
Next, let RSa be the matrix that has as columns the vectors l 1−ρ
(P̃l )j,Sa , l ∈
1/3

ρ
(P̃1 )j,Sa for j ∈ [m] for
[2, L], j ∈ [m] and A is the matrix that has as columns the vectors 1−ρ
some L = O(poly(n)), similarly as in the proof of Theorem 3.1.
We claim that RSa RS⊤a is τ spectrally bounded by ρF F ⊤ .
 
l 1/3
√
ρ
1
Indeed, for any l > 2, we have k l 1−ρ
P̃l k . ρl/3 kP̃l k . ρl/3 mnp Hence,
RSa RS⊤a  ρ2/3 GG⊤ + ρ4/3 HH ⊤ + ρ2 LL⊤ +
 3ρ

2/3

τ (F F

⊤

⊤

+ σmin (F F )) + ρ

X

8/3

38

ρ2l/3 mnp2

l≥4

mnp2 - ρ2/3 τ (F F ⊤ + σmin (F F ⊤ ))

(D.5)

where the first inequality holds since HH ⊤ , GG⊤ , LL⊤ are τ -spectrally bounded bounded by F and
the second since σmin (F F ⊤ ) & np and τ ≥ 1. Let τ ′ = ρ2/3 τ . Since we are assuming the matrix F
is O(1)-incoherent, we can apply Theorem C.10, and claim the output of Algorithm 4 are matrices
Qa , Qb , Qc which are τ -approximate whitening matrices for A, B, C respectively.
By Lemma D.1, we have again
1
1
. (1 + τ ′ )kAA⊤ k, σmin (Qa ) ≥
& (1 − τ ′ )σmin (AA⊤ )
σmax (Qa ) ≤
1 − τ′
1 + τ′

Then, applying Theorem 5.4, we get that we recover vectors (âi , b̂i , ĉi ) are O(η ′ )-close to
1/3
ρ
(1 − exp(Wi )), for all i ∈ [m]. for
1−ρ


3/2
η ′ . max(kQa k , kQb k , kQc k)1/2 · τ ′ + σ −3/2 kEk{1,2},{3}
√
Recall that τ ′ = ρ2/3 τ , and kQa k ≤ ρ2/3 σmax (F ) . ρ2/3 mnp and kEk{1,2},{3} ≤ ρ4 m(np)3/2 and
σ & ρ2/3 np, we obtain that
!
4 m(np)3/2
√
√
ρ
. ρ4/3 mnpτ 3/2
η ′ . ρ1/3 mnp (τ ρ2/3 )3/2 + 2/3
3/2
(ρ np)


where the last inequality holds since ρ3 m = o(1) = o(τ ).
1/3 (â , b̂ , ĉ ) are O(η ′ /ρ1/3 ) = O(η)-close to 1 −
However, this directly implies that ( 1−ρ
i i i
ρ )

exp(−Wi ), i ∈ [m], which in turn implies (ãi , b̃i , c̃i ) are O(η) close to exp(−Wi ).
Argument for recovering Wi from exp(−Wi ) is then exactly the same as the one in Theorem
3.1.

E

Sample complexity and bias of the PMI estimator

Finally, we consider the issue of sample complexity. The estimator we will use for the PMI matrix
will simply be the plug-in estimator, namely:
ˆ i,j = log P̂r[si = 0 ∧ sj = 0]
PMI
ˆ j = 0]
P̂r[si = 0]Pr[s

(E.1)

Notice that this estimator is biased, but as the number of samples grows, the bias tends to zero.
Formally, we can show:
Lemma E.1. If the number of samples N satisfies
1
1
log m
N≥
mini6=j {Pr[si = 0 ∧ sj = 0} δ2

ˆ i,j − PMIi,j | ≤ δ, ∀i 6= j.
with high probability |PMI

Proof. Denoting ∆i,j = P̂r[si = 0 ∧ sj = 0] − Pr[si = 0 ∧ sj = 0] and ∆i = P̂r[si = 0] − Pr[si = 0],
we get that
Pr[si = 0 ∧ sj = 0] + ∆i,j
ˆ i,j = log P̂r[si = 0 ∧ sj = 0] = log
PMI
ˆ
(Pr[s
P̂r[si = 0]Pr[sj = 0]
i = 0] + ∆i )(Pr[sj = 0] + ∆j )






∆i,j
∆i
∆j
= PMIi,j + log 1 +
− log 1 +
− log 1 +
Pr[si = 0 ∧ sj = 0]
Pr[si = 0]
Pr[sj = 0]

39

Furthermore, we have that
2
3x

2x
2+x

≤ log(1 + x) ≤

√x ,
x+1

for x ≥ 0, which implies that when x ≤ 1,

≤ log(1 + x) ≤ x. From this it follows that if


∆i
∆i,j
, max
≤δ
max max
i Pr[si = 0]
i,j Pr[si = 0 ∧ sj = 0]

we have

δ
ˆ i,j ≤ PMIi,j + δ
≤ PMI
3
Note that it suffices to show that if N > 1−4pmax1 mρmax δ12 log m, we have


∆i
∆i
> (1 + δ) ∨
< (1 − δ) ≤ exp(− log2 m)
Pr
Pr[si = 0]
Pr[si = 0]
PMIi,j −

and
Pr




∆i,j
∆i,j
> (1 + δ) ∨
< (1 − δ) ≤ exp(− log2 m)
Pr[si = 0 ∧ sj = 0]
Pr[si = 0 ∧ sj = 0]

since this implies

(E.2)

(E.3)




∆i,j
∆i
max max
, max
≤δ
i,j Pr[si = 0 ∧ sj = 0]
i Pr[si = 0]
with high probability by a simple union bound.
Both (E.2) and (E.3) will follow by a Chernoff bound.
Indeed, consider (E.2) first. We have by Chernoff
s
!
#
"
log N
Pr[si = 0] ≤ exp(− log2 N )
Pr ∆i > 1 +
N Pr[si = 0]
Hence, if N >
2

1
1
Pr[si =0] δ2

log m, we get that 1 − δ ≤

∆i
Pr[si =0]

≤ 1 + δ with probability at least

1 − exp(log m).
The proof of (E.3) is analogous – the only difference being that the requirement is that N >
1
1
Pr[si =0] δ2 log m which gives the statement of the lemma.
Virtually the same proof as above shows that:
Lemma E.2. If the number of samples N satisfies
N≥

1
1
log m
mini6=j6=k {Pr[si = 0 ∧ sj = 0 ∧ sk = 0} δ2

ˆ i,j,k − PMITi,j,k | ≤ δ, ∀i 6= j 6= k.
with high probability |PMIT
As an immediate corollary, we get:
Corollary E.3. If the number of samples N satisfies
N≥

1
1
log m
mini6=j6=k {Pr[si = 0 ∧ sj = 0 ∧ sk = 0} δ2

1
1
log m
mini6=j6=k {Pr[si = 0 ∧ sj = 0 ∧ sk = 0} δ2
with high probability for any equipartition Sa , Sb , Sc
N≥

ˆ Sa ,S ,Sc − PMITSa ,S ,Sc k{1,2},{3} . n3 δ
kPMIT
b
b
40

F

Matrix Perturbation Toolbox

In this section we discuss standard matrix perturbation inequalities. Many results in this section
b = A + E, the perturbation in individual singular
can be found in Stewart and Sun [Ste77]. Given A
values can be bounded by Weyl’s theorem:
b = A+E, we know σk (A)−kEk ≤ σk (A)
b ≤ σk (A)+kEk.
Theorem F.1 (Weyl’s theorem). Given A
For singular vectors, the perturbation is bounded by Wedin’s Theorem:

Lemma F.2 (Wedin’s theorem; Theorem 4.1, p.260 in [Ste90]). Given matrices A, E ∈ Rm×n with
m ≥ n. Let A have the singular value decomposition


Σ1 0
A = [U1 , U2 , U3 ]  0 Σ2  [V1 , V2 ]⊤ .
0
0

b = A + E, with analogous singular value decomposition. Let Φ be the matrix of canonical
Let A
b1 , and Θ be the matrix of canonical angles
angles between the column span of U1 and that of U
between the column span of V1 and that of Vb1 . Suppose that there exists a δ such that
min |[Σ1 ]i,i − [Σ2 ]j,j | > δ,
i,j

and

min |[Σ1 ]i,i | > δ,
i,i

then
k sin(Φ)k2 + k sin(Θ)k2 ≤ 2

kEk2
.
δ2

Perturbation bound for pseudo-inverse When we have a lowerbound on σmin (A), it is easy
to get bounds for the perturbation of pseudoinverse.
Theorem F.3 (Theorem 3.4 in [Ste77]). Consider the perturbation of a matrix A ∈ Rm×n : B =
A + E. Assume that rank(A) = rank(B) = n, then
√
kB † − A† k ≤ 2kA† kkB † kkEk.
Note that this theorem is not strong enough when the perturbation is only known to be τ spectrally bounded in our definition.

41

