On Learning Mixtures of Well-Separated Gaussians

arXiv:1710.11592v1 [] 31 Oct 2017

Oded Regev∗

Aravindan Vijayaraghavan†

Abstract
We consider the problem of efficiently learning mixtures of a large number of spherical
Gaussians, when the components of the mixture are well separated. In the most basic form of
this problem, we are given samples from a uniform mixture of k standard spherical Gaussians
with means µ1 , . . . , µk ∈ Rd , and the goal is to estimate the means up to accuracy δ using
poly(k, d, 1/δ) samples.
In this work, we study the following question: what is the minimum separation needed
between the means for solving this task? The best known algorithm due to Vempala and
1/4
Wang [JCSS 2004] requires a separation of roughly min{k, d} . On the other hand, Moitra
and Valiant [FOCS 2010] showed that with separation o(1), exponentially many samples are
required. We address the significant gap between these two bounds, by showing the following
results.
√
• We show that with separation o( log k), super-polynomially many samples are required.
In fact, this holds even when the k means of the Gaussians are picked at random in
d = O(log k) dimensions.
√
• We show that with separation Ω( log k), poly(k, d, 1/δ) samples suffice. Notice that the
bound on the separation is independent of δ. This result is based on a new and efficient
“accuracy boosting” algorithm that takes as input coarse estimates of the true means and
in time (and samples) poly(k, d, 1/δ) outputs estimates of the means√up to arbitrarily
good
√
accuracy δ assuming the separation between the means is Ω(min{ log k, d}) (independently of δ). The idea of the algorithm is to iteratively solve a “diagonally dominant”
system of non-linear equations.
We√also (1) present a computationally efficient algorithm in d = O(1) dimensions with only
Ω( d) separation, and (2) extend our results to the case that components might have different
weights and variances. These results together essentially characterize the optimal order of
separation between components that is needed to learn a mixture of k spherical Gaussians with
polynomial samples.

1

Introduction

Gaussian mixture models are one of the most widely used statistical models for clustering. In this
model, we are given random samples, where each sample point x ∈ Rd is drawn independently from
∗

Courant Institute of Mathematical Sciences, New York University. Supported by the Simons Collaboration on
Algorithms and Geometry and by the National Science Foundation (NSF) under Grant No. CCF-1320188.
†
Department of Electrical Engineering and Computer Science, Northwestern University. Supported by the National
Science Foundation (NSF) under Grant No. CCF-1652491 and CCF-1637585. Any opinions, findings, and conclusions
or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of
the NSF.

1

one of k Gaussian components according to mixing weights w1 , w2 , . . . , wk , where each Gaussian
component j ∈ [k] has a mean µj ∈ Rd and a covariance Σj ∈ Rd×d . We focus on an important
special case of the problem where each of the components is a spherical Gaussian, i.e., the covariance
matrix of each component is a multiple of the identity. If f represents the p.d.f. of the Gaussian
mixture G, and gj represents the p.d.f. of the jth Gaussian component,
k

X

1
wj gj (x).
gj = d exp −πkx − µj k22 /σj2 , f (x) =
σj
j=1
The goal is to estimate the parameters {(wj , µj , σj ) : j ∈ [k]} up to required accuracy δ > 0 in time
and number of samples that is polynomial in k, d, 1/δ.
Learning mixtures of Gaussians has a long and rich history, starting with the work of Pearson
[38]. (See Section 1.2 for an overview of prior work.) Most of the work on this problem, especially
in the early years but also recently, is under the assumption that there is some minimum separation
between the means of the components in the mixture. Starting with work by Dasgupta [16], and
continuing with a long line of work (including [4, 45, 2, 29, 39, 17, 14, 30, 7, 8, 46, 19]), efficient
algorithms were found under mild separation assumptions. Considering for simplicity the case of
uniform mixtures (i.e., all weights are 1/k) of standard Gaussians (i.e., spherical with σ = 1), the
best known result due to Vempala and Wang [45] provides an efficient algorithm (both in terms
of samples and running time) under separation of at least min{k, d}1/4 poly log(dk/δ) between any
two means.
A big open question in the area is whether efficient algorithms exist under weaker separation
assumptions. It is known that when the separation is o(1), a super-polynomial number of samples
is required (e.g., [35, 3, 26]), but the gap between this lower bound and the above upper bound of
min{k, d}1/4 poly log(dk/δ) is quite wide. Can it be that efficient algorithms exist under only Ω(1)
separation? In fact, prior to this work, this was open even in the case of d = 1.
Question 1.1. What is the minimum order of separation that is needed to learn the parameters
of a mixture of k spherical Gaussians up to accuracy δ using poly(d, k, 1/δ) samples?

1.1

Our Results

By improving both the lower bounds and the upper bounds mentioned above, we characterize
(up to constants) the minimum separation needed to learn the mixture from polynomially many
samples.
Our first result shows super-polynomial lower bounds when the separation is of the order
√
o( log k). In what follows, ∆param (G, G̃) represents the “distance” between the parameters of the
two mixtures of Gaussians G, G̃ (see Definition 2.2 for the precise definition).
√
Informal Theorem 1.2 (Lower Bounds). For any γ(k) = o( log k), there are two uniform mixtures of standard spherical Gaussians G, G̃ in d = O(log k) dimensions with means {µ1 , . . . , µk }, {µ̃1 , µ̃2 , . . . , µ̃k }
respectively, that are well separated
∀i 6= j ∈ [k] : kµi − µj k2 ≥ γ(k), and kµ̃i − µ̃j k2 ≥ γ(k),
and whose parameter distance is large ∆param ({µ1 , . . . , µk }, {µ̃1 , . . . , µ̃k }) = Ω(1), but have very
small statistical distance kG − G̃kT V ≤ k−ω(1) .
2

The above statement implies that we need at least kω(1) many samples to distinguish between
G, G̃, and identify G. See Theorem 3.1 for a formal statement of the result. In fact, these sample
complexity
√ lower bounds hold even when the means of the Gaussians are picked randomly in a ball
of radius d in d = o(log k) dimensions. This rules out obtaining smoothed analysis guarantees for
small dimensions (as opposed to [10, 3] which give polytime algorithms for smoothed mixtures of
Gaussians in kΩ(1) dimensions).
√
Our next result shows that the separation of Ω( log k) is tight – this separation suffices to learn
the parameters of the mixture with polynomial samples. We state the theorem for the special case
of uniform mixtures of spherical Gaussians. (See Theorem 5.1 for the formal statement.)
Informal Theorem 1.3 (Tight Upper Bound in terms of k). There exists a universal constant
c > 0, such that given samples from a uniform mixture of standard spherical Gaussians in Rd with
well-separated means, i.e.,
p
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ c log k
(1)
there is an algorithm that for any δ > 0 uses only poly(k, d, 1/δ) samples and with high probability
finds {µ̃1 , µ̃2 , . . . , µ̃k } satisfying ∆param ({µ1 , . . . , µk }, {µ̃1 , . . . , µ̃k }) ≤ δ.

While the above algorithm uses only poly(k, d, 1/δ) samples, it is computationally inefficient.
Our next result shows that in constant dimensions, one can obtain a computationally efficient
algorithm. In fact, in such low dimensions a separation of order Ω(1) suffices.
Informal Theorem 1.4 (Efficient algorithm in low dimensions). There exists a universal constant
c > 0, such that given samples from a uniform mixture of standard spherical Gaussians in Rd with
well-separated means, i.e.,
√
(2)
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ c d
there is an algorithm that for any δ > 0 uses only polyd (k, 1/δ) time (and samples) and with high
probability finds {µ̃1 , µ̃2 , . . . , µ̃k } satisfying ∆param ({µ1 , . . . , µk }, {µ̃1 , . . . , µ̃k }) ≤ δ.
See Theorem 6.1 for a formal statement. An important feature of the above two algorithmic
results is that the separation is independent of the accuracy δ that we desire in parameter estimation
(δ can be arbitrarily small compared to k and d).
These results together almost give a tight
characterization (up to constants) for the amount of separation needed to learn with poly(k, d, 1/δ)
samples.
Iterative Algorithm. The core technical portion of Theorem 1.3 and Theorem 1.4 is a new
iterative algorithm, which is the main algorithmic contribution of the paper. This algorithm takes
coarse estimates of the means, and iteratively refines them to get arbitrarily good accuracy δ. We
now present an informal statement of the guarantees of the iterative algorithm.
Informal Theorem 1.5 (Iterative Algorithm Guarantees). There exists a universal constant c > 0,
such that given samples from a uniform mixture of standard spherical Gaussians in Rd with wellseparated means, i.e.
√
p
(3)
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ c min{ log k, d}
3

and suppose we are given initializers µ̃1 , . . . , µ̃k for the means µ1 , . . . , µk satisfying
∀j ∈ [k],


1
kµj − µ̃j k2 ≤ 1/poly min{d, k} .
σj

There exists an iterative algorithm that for any δ > 0 that runs in poly(k, d, 1/δ) time (and sam(T )
(T )
ples), and after T = O(log log(k/δ)) iterations, finds with high probability µ1 , . . . , µk such that
(T )
(T )
∆param ({µ1 , . . . , µk }, {µ1 , . . . , µk }) ≤ δ.
The above theorem also holds when the weights and variances are unequal. See Theorem 4.1
for a formal statement. Note that in the above result, the desired accuracy δ can be arbitrarily
small compared to k, and the separation required does not depend on δ. To prove the polynomial
identifiability results (Theorems 1.3 and 1.4), we first find coarse estimates of the means that
serve as initializers to this iterative algorithm, which then recovers the means up to arbitrarily fine
accuracy independent of the separation.
The algorithm works by solving a system of non-linear equations that is obtained by estimating
simple statistics (e.g., means) of the distribution restricted to certain carefully chosen regions. We
prove that the system of non-linear equations satisfies a notion of “diagonal dominance” that allows
us to leverage iterative algorithms like Newton’s method and achieve rapid (quadratic) convergence.
The techniques developed here can find such initializers using only poly(k, d) many samples,
but use time that is exponential in k. This leads to the following natural open question:
Open Question 1.6. Given a mixture of spherical Gaussians with equal weights and variances,
and with separation
p
∀i 6= j ∈ [k], kµi − µj k2 ≥ c log k

for some sufficiently large absolute constant c > 0, is there an algorithm that recovers the parameters
up to δ accuracy in time poly(k, d, 1/δ)?
Our iterative algorithm shows that to resolve this open question affirmatively, it is enough to
find initializers that are reasonably close to√the true parameters. In fact, a simple amplification
argument shows that initializers that are c log k/8 close to the true means will suffice for this
approach.
Our iterative algorithm is reminiscent of some commonly used iterative heuristics, such as
Lloyd’s Algorithm and especially Expectation Maximization (EM). While these iterative methods
are the practitioners’ method-of-choice for learning probabilistic models, they have been notoriously
hard to analyze. We believe that the techniques developed here may also be useful to prove
guarantees for these heuristics.

1.2

Prior Work and Comparison of Results

Gaussian mixture models are among the most widely used probabilistic models in statistical inference [38, 41, 42]. Algorithmic results fall into two broad classes — separation-based results, and
moment-based methods that do not assume explicit geometric separation.
Separation-based results. The body of work that is most relevant to this paper assumes that
there is some minimum separation between the means of the components in the mixture. The first
polynomial time algorithmic guarantees for mixtures of Gaussians were given by Dasgupta [16],
4

who showed how to learn mixtures of spherical Gaussians when the separation is of the order of
d1/2 . This was later improved by a series of works [4, 45, 2, 29, 17, 14, 30, 7] for both spherical
Gaussians and general Gaussians. The algorithm of Vempala and Wang [45] gives the best known
result, and uses PCA along with distance-based clustering to learn mixtures of spherical Gaussians
with separation
kµi − µj k2 ≥ (min{k, d}1/4 log1/4 (dk/δ) + log1/2 (dk/δ))(σi + σj ).
We note that all these clustering-based algorithms require a separation that either implicitly or explicitly depend on the estimation accuracy δ.1 Finally, although not directly related to our work, we
note that a similar separation condition was shown to suffice also for non-spherical Gaussians [14],
where separation is measured based on the variance along the direction of the line joining the respective means (as opposed, e.g., to the sum of maximum variances kΣi k + kΣj k which could be
much larger).
Iterative methods like Expectation Maximization (EM) and Lloyd’s algorithm (sometimes called
the k-means heuristic) are commonly used in practice to learn mixtures of spherical Gaussians but,
as mentioned above, are notoriously hard to analyze. Dasgupta and Schulman [17] proved that a
variant of the EM algorithm learns mixtures of spherical Gaussians with separation of the order
of d1/4 polylog(dk). Kumar and Kannan [30] and subsequent work [7] showed that the spectral
clustering heuristic (i.e., PCA followed by Lloyd’s algorithm) provably recovers the clusters in a
rather wide family of distributions which includes non-spherical
√ Gaussians; in the special case of
spherical Gaussians, their analysis requires separation of order k.
Very recently, the EM algorithm was shown to succeed for mixtures of k = 2 spherical Gaussians
with Ω(σ) separation [8, 46, 19] (we note that in this setting with k = O(1), polynomial time
guarantees are also known using other algorithms like the method-of-moments [28], as we will see
in the next paragraph). SDP-based algorithms have also been studied in the context of learning
mixtures of spherical Gaussians with a similar separation requirement of Ω(k maxi σi ) [34]. The
question of how much separation between the components is necessary was also studied empirically
by Srebro et al. [39], who observed that iterative heuristics successfully learn the parameters under
much smaller separation compared to known theoretical bounds.
Moment-based methods. In a series of influential results, algorithms based on the method-ofmoments were developed by [28, 35, 9] for efficiently learning mixtures of k = O(1) Gaussians under
arbitrarily small separation. To perform parameter estimation up to accuracy δ, the running time
2
of the algorithms is poly(d, 1/wmin , 1/δ)O(k ) (this holds for mixtures of general Gaussians). This
exponential dependence on k is necessary in general, due to statistical lower bound results [35].
The running time dependence on δ was improved in the case of k = 2 Gaussians in [26].
Recent work [27, 11, 25, 10, 3, 24] use uniqueness of tensor decompositions (of order 3 and
above) to implement the method of moments and give polynomial time algorithms assuming the
means are sufficiently high dimensional, and do not lie in certain degenerate configurations. Hsu and
Kakade [27] gave a polynomial time algorithm based on tensor decompositions to learn a mixture of
spherical Gaussians, when the means are linearly independent. This was extended by [25, 10, 3] to
give smoothed analysis guarantees to learn “most” mixtures of spherical Gaussians when the means
1
Such a dependency on δ seems necessary for clustering-based algorithms that cluster every point accurately with
high probability.

5

are in d = kΩ(1) dimensions. These algorithms do not assume any strict geometric separation conditions and learn the parameters in poly(k, d, 1/δ) time (and samples), when these non-degeneracy
assumptions hold. However, there are many settings where the Gaussian mixture consists of many
clusters in a low dimensional space, or have their means lying in a low-dimensional subspace or
manifold, where these tensor decomposition guarantees do not apply. Besides, these algorithms
based on tensor decompositions seem less robust to noise than clustering-based approaches and
iterative algorithms, giving further impetus to the study of the latter algorithms as we do in this
paper.
Lower Bounds. Moitra and Valiant [35] showed that exp(k) samples are needed to learn the
parameters of a mixture of k Gaussians√[35]. In fact, the lower bound instance of [35] is one
dimensional, with separation of order 1/ k. Anderson et al. [3] proved a lower bound on sample
complexity that is reminiscent of our Theorem 1.2. Specifically, they obtain a super-polynomial
lower bound assuming separation O(σ/poly log(k)) for d = O(log k/
√ log log k). This is in contrast
to our lower bound which allows separation greater than σ, or o(σ log k) to be precise.
Other Related work. While most of the previous work deals with estimating the parameters
(e.g., means) of the Gaussians components in the given mixture G, a recent line of work [23, 40, 18]
focuses on the task of learning a mixture of Gaussians G ′ (with possibly very different parameters)
that is close in statistical distance i.e., kG − G ′ kT V < δ (this is called “proper learning”, since the
hypothesis that is output is also a mixture of k Gaussians). When identifiability using polynomial
samples is known for the family of distributions, proper learning automatically implies parameter
estimation. Algorithms for proper learning mixtures of spherical Gaussians [40, 18] give polynomial
sample complexity bounds when d = 1 (note that the lower bounds of [35] do not apply here) and
have run time dependence that is exponential in k; the result of [23] has sample complexity that
is polynomial in d but exponential in k. Algorithms that take time and samples poly(k, 1/δ)d
are also known for “improper learning” and density estimation for mixtures of k Gaussians (the
hypothesis class that is output may not be a mixture of k Gaussians) [15, 12]. We note that known
algorithms have sample complexity that is either exponential in d or k, even though proper learning
and improper learning are easier tasks than parameter estimation. To the best of our knowledge,
better bounds are not known under additional separation assumptions.
A related problem in the context of clustering graphs and detecting communities is the problem
of learning a stochastic block model or planted partitioning model [33]. Here, a sharp threshold
phenomenon involving an analogous separation condition (between intra-cluster probability of edges
and inter-cluster edge probability) is known under various settings [36, 32, 37] (see the recent survey
by Abbe [1] for details). In fact, the algorithm of Kumar and Kannan [30] give a general separation
condition that specializes to separation between the means for mixtures of Gaussians, and separation
between the intra-cluster and inter-cluster edge probabilities for the stochastic block model.

1.3

Overview of Techniques

√
Lower bound for O( log k) separation. The sample complexity lower bound proceeds by
showing a more general statement: in any large enough collection of uniform mixtures, for all but
a small fraction of the mixtures, there is at least one other mixture in the collection that is close in
statistical distance (see Theorem 3.2). For our lower bounds, we will just produce a large collection
of uniform mixtures of well-separated spherical Gaussians in d = c log k dimensions, whose pairwise
6

parameter distances are reasonably
large. In fact, we can even pick the means of these mixtures
√
randomly in a ball of radius d in d = c log k dimensions; w.h.p. most of these mixtures will need
at least kω(1) samples to identify.
To show the above pigeonhole style statement about large collections of mixtures, we will
associate with a uniform mixture having means µ1 , . . . , µk , the following quantities that we call
“mean moments,” and we will use them as a proxy for the actual moments of the distribution:
(M1 , . . . , MR ) where ∀1 ≤ r ≤ R : Mr =

k
1 X ⊗r
µj .
k
j=1

The mean moments just correspond to the usual moments of a mixture of delta functions centered
at µ1 , . . . , µk . Closeness in the first R = O(1/ε) mean moments (measured in injective tensor norm)
implies that the two corresponding distributions are ε close in statistical distance (see Lemma 3.7
and Lemma 3.8). The key step in the proof uses a careful packing argument to show that for
most mixtures in a large enough collection, there is a different mixture in the collection that
approximately matches in the first R mean moments (see Lemma 3.6).
Iterative Algorithm. Our iterative algorithm
will function in both settings of interest: the
√
log
k)
separation,
and the low-dimensional setting when
high-dimensional setting when
we
have
Ω(
√
d < log k and we have Ω( d) separation. For the purpose of this description, let us assume δ is
arbitrarily small compared to d and k. In our proposed algorithm, we will consider distributions
obtained by restricting the support to just certain regions around the initializers z1 = µ̃1 , . . . , zk =
µ̃k that are somewhat close to the means µ1 , µ2 , . . . , µk respectively. Roughly speaking, we first
partition the space into a Voronoi partition given by {zj : j ∈ [k]}, and then for each component
j ∈ [k] in G, let Sj denote the region containing zj (see Definition 4.3 for details). For each j ∈ [k]
we consider only the samples in the set Sj and let uj ∈ Rd be the (sample) mean of these points in
Sj , after subtracting zj .
The regions are chosen in such a way that Sj has a large fraction of the probability mass from
the jth component, and the total probability
mass from the other components is relatively small (it
√
will be at most 1/poly(k) with Ω( log k) separation, and Od (1) with Ω(1) separation in constant
dimensions). However, since δ can be arbitrarily small functions of k, d, there can still be a relatively
large contribution from the other components. For instance, in the low-dimensional case with O(1)
separation, there can be Ω(1) mass from a single neighboring component! Hence, uj does
p not give
a δ-close estimate for µj (even up to scaling),
√ unless the separation is at least of order log(1/δ)
– this is too large when δ = k−ω(1) with log k separation, or δ = od (1) with Ω(1) separation in
constant dimensions.
Instead we will use these statistics to set up a system of non-linear equations where the unknowns
are the true parameters and solve for them using the Newton method. We will use the initializers
(0)
zj = µj , to define the statistics that give our equations. Hence the unknown parameters {µi : i ∈
[k]} satisfy the following equation for each j ∈ [k]:
k
X
i=1

wi

Z

y∈Sj

(y

− zj ) · σj−d exp



πky − µi k22
−
σi2



dy = uj .

(4)

Note that in the above equation, the only unknowns or variables are the true means {µi : i ∈ [k]}.
After scaling the equations, and a suitable change of variables xj = µj /σj to make the system
7

“dimensionless” we get a non-linear system of equations denoted by F (x) = b. For the above
system, x∗i = µi /σi represents a solution to the system given by the parameters of G. The Newton
algorithm uses the iterative update
x(t+1) = x(t) + (F ′ (x(t) ))−1 (b − F (x(t) )).
For the Newton method we need access to the estimates for b, and the derivative matrix F ′ (the
Jacobian) evaluated at x(t) . The derivative of the j equation w.r.t. xi corresponds to
Z
2πwi
(y − zj )(y − σi xi )T gσi xi ,σi (y) dy ,
∇xi Fj (x) =
wj σj σi y∈Sj
where gσi xi ,σi (y) represents the p.d.f. at a point y due to a spherical Gaussian with mean at σi xi
and covariance σi2 /(2π) in each direction. Unlike usual applications of the Newton method, we do
not have closed form expressions for F ′ (the Jacobian), due to our definition of the set Sj . However,
we will instead be able to estimate the Jacobian at x(t) by calculating the above expression (RHS)
(t)
by considering a Gaussian with mean σi xi and variance σi2 /(2π). The Newton method can be
shown to be robust to errors in b, F, F ′ (see Theorem B.2).
We want to learn each of the k means up to good accuracy; hence we will measure the error and
convergence in k·k∞ norm. This
√ is important in low dimensions since measuring convergence in
ℓ2 norm will introduce extra k factors, that are prohibitive for us since the means are separated
only by Θd (1). The convergence of the Newton’s method depends on upper bounding the operator
norm of the inverse of the Jacobian k(F ′ )−1 k and the second-derivative kF ′′ k, with the initializer
being chosen δ-close to the true parameters so that δk(F ′ )−1 kkF ′′ k < 1/2.
The main technical effort for proving convergence is in showing that the inverse (F ′ )−1 evaluated
at any point in the neighborhood around x∗ is well-conditioned. We will show the convergence of
the Newton method by showing “diagonal dominance” properties of the dk × dk matrix F ′ . This
uses the separation between
√ the means of the components, and the properties of the region Sj that
we have defined. For Ω( log k) separation, this uses standard facts about Gaussian concentration
to argue that each of the (k − 1) off-diagonal blocks (in the jth row of F ′ ) is at most 1/(2k) factor of
the corresponding diagonal term. With Ω(1) separation in d = O(1) dimensions, we can not hope
to get such a uniform bound on all the off-diagonal blocks (a single off-diagonal block can itself be
Ωd (1) times the corresponding diagonal entry). We will instead use careful packing arguments to
show that the required diagonal dominance condition (see Lemma 4.13 for a statement).
Hence, the initializers are used to both define the regions Sj , and as initialization for the
Newton method. Using this diagonal dominance in conjunction with initializers (Theorem 5.2) and
the (robust) guarantees of the Newton method (Theorem B.2) gives rapid convergence to the true
parameters.

2

Preliminaries

Consider a mixture of k spherical Gaussians G in Rd that has parameters {(wj , µj , σj ) : j ∈ [k]}. The
jth component has mean µj and covariance σj2 /2π · Id×d . For µ ∈ Rd , σ ∈ R+ , let gµ,σ : Rd → R+
represent the p.d.f. of a spherical Gaussian centered at µ and with covariance σ 2 /(2π) · Id×d . We
will use f to represent the p.d.f. of the mixture of Gaussians G, and gj to represent the p.d.f. of
the jth Gaussian component.
8

Definition 2.1 (Standard mixtures of Gaussians). A standard mixture of k Gaussians with means
µ1 , . . . , µk ∈ Rd is a mixture of k spherical Gaussians {( k1 , µj , 1) : j ∈ [k]}.
A standard mixture is just a uniform mixture of spherical Gaussians with all covariances σ 2 =
1/(2π). Before we proceed, we define the following notion of parameter “distance” between mixtures
of Gaussians:
Definition 2.2 (Parameter distance). Given two mixtures of Gaussians in Rd , G = {(wj , µj , σj ) :
j ∈ [k]} and G ′ = {(wj′ , µ′j , σj′ ) : j ∈ [k]}, define
∆param G, G

′



′
k kµ − µ′
k
k
X
X
X
|σj − σπ(j)
|
|wj − wπ(j) |
j
π(j) k2
= min
+
+
.
′
′
π∈Permk
min{wj , wπ(j) }
min{σj , σπ(j) }
min{σj , σπ(j) }
j=1

j=1

j=1

For standard mixtures, the definition simplifies to

∆param (µ1 , . . . , µk ), (µ′1 , . . . , µ′k ) =

min

π∈Permk

k
X
kµj − µ′π(j) k2 .
j=1

Note that this definition is invariant to scaling the variances (for convenience). We note that
parameter distance is not a metric, but it is just a convenient way of measure closeness of parameters
between two distributions.
The distance between two individual Gaussian components can also be measured in terms of
the total variation distance between the components
[35]. For instance, in the case of standard
√
spherical Gaussians, a parameter distance of c log k corresponds to a total variation distance of
2
k−O(c ) .
Definition 2.3 (ρ-bounded mixtures). For ρ ≥ 1, a mixture of spherical Gaussians G = {(wj , µj , σj )}kj=1
in Rd is called ρ-bounded if for each j ∈ [k], kµj k2 ≤ ρ and ρ1 ≤ σj ≤ ρ. In particular, a standard
mixture is ρ-bounded if for each j ∈ [k], kµj k2 ≤ ρ.
Also, for a given mixture of k spherical gaussians G = {(wj , µj , σj ) : j ∈ [k]}, we will denote
wmin = minj∈[k] wj , σmax = maxj∈[k] σj and σmin = minj∈[k] σj .
In the above notation the bound ρ can be thought of as a sufficiently large polynomial in k,
since we are aiming for bounds that are polynomial in k. Since we can always scale the points by
an arbitrary factor without affecting the performance of the algorithm, we can think of ρ as the
(multiplicative) range of values taken by the parameters {µi , σi : i ∈ [k]}. Since we want separation
bounds independent of k, we will denote individual aspect ratios for variances and weights given
by ρσ = maxi∈[k] σi / mini∈[k] σi , and ρw = maxi∈[k] wi / mini∈[k] wi .
Finally, we list some of the conventions used in this paper. We will denote by N (0, σ 2 ) a normal
random variable with mean 0 and variance σ 2 . For x ∈ R generated according to N (0, σ 2 ), let Φ̃0,σ (t)
denote the probability that x > t, and let Φ̃−1
0,σ (y) denote the quantile t at which Φ̃0,σ (t) ≤ y. For
d
′
any function f : R → R, f will denote the first derivative (or gradient)
of the function, and f ′′
R
will denote the second derivative (or Hessian). We define kf k1,S = S |f (x)|dx to be the L1 norm of
f restricted to the set S. Typically, we will use indices i, j to represent one of the k components of
the mixture, and we will use r (and s) for coordinates. For a vector x ∈ Rd , we will use x(r) denote
the rth coordinate. Finally, we will use w.h.p. in statements about the success of algorithms to
represent probability at least 1 − γ where γ = (d + k)−Ω(1) .
9

Norms.

For any p ≥ 1, given a matrix M ∈ Rd×d , we define the matrix norm
kM kp→p =

2.1

max

x∈Rd :kxkp =1

kM xkp .

Notation and Preliminaries about Newton’s method

Consider a system of m non-linear equations in variables u1 , u2 , . . . , um :
∀j ∈ [m], fj (u1 , . . . , um ) = bj .
Let F ′ = J(u) ∈ Rm×m be the Jacobian of the system given by the non-linear functional
∂fj (u)
f : Rm → Rm , where the (j, i)th entry of J is the partial derivative ∂u
is evaluated at u.
i
(0)
Newton’s method starts with an initial point u , and updates the solution using the iteration:

−1 

u(t+1) = u(t) + J(u(t) )
bj − f (u(t) ) .

Standard results shows quadratic convergence of the Newton method for general normed spaces [5].
We restrict our attention in the restricted setting where both the range and domain of f is Rm ,
equipped with an appropriate norm k·k to measure convergence.
Theorem 2.4 (Theorem 5.4.1 in [5]). Assume u∗ ∈ Rm is a solution to the equation f (y) = b
where f : Rm → Rm and the inverse Jacobian J −1 exists in a neighborhood N = {u : ku − u∗ k ≤
ku(0) − u∗ k}, and F ′ : Rm → Rm×m is locally L-Lipschitz continuous in the neighborhood N i.e.,
∀u, v ∈ N, kF ′ (u) − F ′ (v)k ≤ Lku − vk. Then we have ku(t+1) − u∗ k ≤ L · kJ(u(t) )−1 k · ku(t) − u∗ k2 .
In particular, for Newton’s method to work, ku0 − u∗ k ≤ (L maxu∈N kJ(u)−1 k)−1 will guarantee convergence. A statement of the robust convergence of Newton’s method in the presence of
estimates is given in Theorem B.2 and Corollary B.4.
We want to learn each of the k sets of parameters up to good accuracy; hence we will measure
the error in ℓ∞ norm. To upper bound kJ −1 k∞→∞ , we will use diagonal dominance properties of
the matrix J. Note that kAk∞→∞ is just the maximum ℓ1 norm of the rows of A. The following
lemma bound kA−1 k∞→∞ for a diagonally dominant matrix A.
Lemma 2.5 ([44]). Consider any square matrix A of size n × n satisfying
X
∀i ∈ [n] aii −
|aij | ≥ α.
j6=i

Then, kA−1 k∞→∞ ≤ 1/α.

√
Lower Bounds with O( log k) Separation

3

Here we show a sample complexity lower bound for learning
standard mixtures of k spherical
√
Gaussians even when the separation is of the order of log k. In fact, this lower bound will also
hold for a random mixture of Gaussians in d ≤ c · log k dimensions (for sufficiently small constant
c) with high probability.2
2

In particular, this rules out polynomial-time smoothed analysis guarantees of the kind shown for d = kΩ(1) in
[10, 3].

10

Theorem 3.1. For any large enough C there exist c, c2 > 0, such that the following holds for
all k ≥ C 8 . Let D be the distribution over standard mixtures of k spherical Gaussians
obtained
√
by picking each of the k means independently and uniformly from a ball of radius d around the
origin in d = c log k dimensions. Let {µ1 , µ2 , . . . , µk } be a mixture chosen according to D. Then
with probability at least 1 − 2/k there exists another √
standard mixture of k spherical Gaussians with
means {µ̃1 , µ̃2 , . . . , µ̃k } such that both mixtures are d bounded and well separated, i.e.,
p
p
∀i, j ∈ [k], i 6= j : kµi − µj k ≥ c2 log k and kµ̃i − µ̃j k ≥ c2 log k,

and their p.d.f.s satisfy

kf − f˜k1 ≤ k−C
(5)
√
even though their parameter distance is at least c2 log k. Moreover, we can take c = 1/(4 log C)
and c2 = C −24 .

Remark. In Theorem 3.1, there is a trade-off between getting a smaller statistical distance ε = k−C ,
and a larger separation between the means in the Gaussian mixture. √When C = ω(1), with
c1 , c = o(1) we see that kf − f˜k1 ≤ k−ω(1) when the separation is o( log k)σ. On the other
′
hand, we can also set C = kε (for some small constant ε′ > 0) to get lower bounds for mixtures
of spherical Gaussians in d = 1 dimension with kf − f˜k1 = exp(−kΩ(1) ) and separation 1/k O(1)
between the means. We note that in this setting of parameters, our lower bound is nearly identical
to that in [35]. Namely, they achieve statistical distance kf − f˜k1 = exp(−Ω(k)) (which is better
than ours) with a similar separation of k−O(1) in one dimension. One possible advantage of our
bound is that it holds with a random choice of means, unlike their careful choice of means.

3.1

Proof of Theorem 3.1

The key to the proof of Theorem 3.1 is the following pigeonhole statement, which can be viewed
as a bound on the covering number (or equivalently, the metric entropy) of the set of Gaussian
mixtures.
Theorem 3.2. Suppose we are√given a collection F of standard
mixtures of spherical Gaussians
√
in d dimensions that are ρ = d bounded, i.e., kµj k ≤ d for all j ∈ [k]. There are universal
constants c0 , c1 ≥ 1, such that for any η > 0, ε ≤ exp(−c1 d), if

 
1
log(1/ε) d
|F| > exp c0
· log(1/ε) log(5d) ,
(6)
η
d

then for at least (1 − η) fraction of the mixtures {µ1 , µ2 , . . . , µk } from F, there is another mixture
{µ̃1 , µ̃2 , . . . , µ̃k } from F with p.d.f. f˜ such that kf − f˜k1 ≤ ε. Moreover, we can take c0 = 8πe and
c1 = 36.
Remark 3.3. Notice that k plays no role in the statement above. In fact, the proof also holds for
mixtures with arbitrary number of components and arbitrary weights.
Claim 3.4. Let x1 , . . . , xN be chosen independently and uniformly from the ball of radius r in Rd .
Then for any 0 < γ < 1, with probability at least 1−N 2 γ d , we have that for all i 6= j, kxi −xj k ≥ γr.

Proof. For any fixed i 6= j, the probability that kxi − xj k ≥ γr is at most γ d , because the volume
of a ball of radius γr is γ d times that of a ball of radius r. The claim now follows by a union
bound.
11

Proof of Theorem 3.1. Set γ := 2−6/c , and consider the following probabilistic procedure. We first
√
let X be a set of (1/γ)d/3 points chosen independently and uniformly from the ball of radius d.
We then output a mixture chosen uniformly from the collection F, defined as the collection of all
standard mixtures of spherical Gaussians obtained by selecting k distinct means from X . Observe
that the output of this procedure is distributed according to D. Our goal is therefore to prove that
with probability at least 1 − 2/k, the output of the procedure satisfies the property in the theorem.
First, by Claim√3.4, with probability at least 1 − γ d/3 ≥ 1 − 1/k, any two points in X are√at
distance at least γ d. It follows that in this case, the means in any mixture in F are at least
√ γ d
apart, and also that any two distinct mixtures in F have a parameter distance of at least γ d since
they must differ in at least one of the means. Note that γ = C −24 for our choice of c, γ.
To complete the proof, we notice that by our choice of parameters, and denoting ε = k−C ,
   
 

|X |
log(1/ε) d
1 dk/3 −k
k
|F| =
≥
· k = k ≥ k · exp c0
· log(1/ε) log(5d) .
k
γ
d
1
The last inequality follows since for our choice of ε = k−C , c = 4 log
C and C is large enough with
C ≥ c0 , so that


√
√
log(1/ε) d
= kc log(C/c) < k, and c0 log(1/ε) log(5d) ≤ c0 C log k log(5c log k) < k.
d

Hence applying Theorem 3.2 to F, for at least 1−1/k fraction of the mixtures in F, there is another
mixture in F that is ε close in total variation distance. We conclude that with probability at least
1 − 2/k, a random mixture in F satisfies all the required properties, as desired.

3.2

Proof of Theorem 3.2

It will be convenient to represent the p.d.f. f (x) of the standard mixture of spherical Gaussians
with means µ1 , µ2 , . . . , µk as a convolution of a standard mean zero Gaussian with a sum of delta
functions centered at µ1 , µ2 , . . . , µk ,


k
X
1
2
δ(x − µj ) ∗ e−πkxk2 .
f (x) = 
k
j=1

Instead of considering the moments of the mixture of Gaussians, we will consider moments of
just the corresponding mixture of delta functions at the means. We will call them “mean moments,”
and we will use them as a proxy for the actual moments of the distribution.
k
1 X ⊗r
µj .
(M1 , . . . , MR ) where ∀1 ≤ r ≤ R : Mr =
k
j=1

To prove Theorem 3.2 we will use three main steps. Lemma 3.6 will show using the pigeonhole
principle that for any large enough collection of Gaussian mixtures F, most Gaussians mixtures
in the family have other mixtures which approximately match in their first R = O(log(1/ε)) mean
moments. This closeness in moments will be measured using the symmetric injective tensor norm
ℓ
defined for an order-ℓ tensor T ∈ Rd as
kT k∗ = max |hT, y ⊗ℓ i|.
y∈Rd
kyk=1

12

Next, Lemma 3.7 shows that the two distributions that are close in the first R mean moments
are also close in the L2 distance. This translates to small statistical distance between the two
distributions using Lemma 3.8.
We will use the following standard packing claim.
Claim 3.5. Let k · k be an arbitrary norm on RD . If x1 , . . . , xN ∈ RD are such that kxi k ≤ ∆ for
all i, and for all i 6= j, kxi − xj k > δ, then N ≤ (1 + 2∆/δ)D . In particular, if x1 , . . . , xN ∈ RD
are such that kxi k ≤ ∆ for all i, then for all but (1 + 2∆/δ)D of the indices i ∈ [N ], there exists a
j 6= i such that kxi − xj k ≤ δ.
Proof. Let K be the unit ball of the norm k · k. Then by assumption, the sets xi + δK/2 for i ∈ [N ]
are disjoint. But since they are all contained in (∆ + δ/2)K,
N≤

Vol((∆ + δ/2)K)
= (1 + 2∆/δ)D .
Vol(δK/2)

The “in particular” part follows by taking a maximal set of δ-separated points.
Lemma 3.6. Suppose we are given√a set F of standard mixtures of spherical Gaussians in d dimen- 
sions with means of length at most d. Then for any integer R ≥ d, if |F| > η1 ·exp (2eR/d)d R log(5d) ,
it holds that for at least (1 − η) fraction of the mixtures {µ1 , µ2 , . . . , µk } in F, there is another mixture {µ̃1 , µ̃2 , . . . , µ̃k } in F satisfying that for r = 1, . . . , R,
k
k
1 X ⊗r 1 X
µj −
(µ̃j )⊗r
k
k
j=1

∗

j=1

≤ (d + 1)−R/4 .

(7)

Proof. With any choice of means µ1 , µ2 , . . . , µk ∈ Rd we can associate a vector of moments
ψ(µ1 , µ2 , . . . , µk ) = (M1 , . . . , MR ) where for r = 1, . . . , R,
k
1 X ⊗r
µj .
Mr =
k
j=1


Notice that the image of ψ lies in a direct sum of symmetric subspaces whose dimension is D = d+R
R
(i.e., the number of ways of distributing R identical balls into (d + 1) different bins). Since R ≥ d,

 



d+R
e(d + R) d
2eR d
D=
≤
≤
.
(8)
R
d
d
We define a norm on these vectors ψ(µ1 , . . . , µk ) in terms of the maximum injective norm of the
mean moments,
kψ(µ1 , . . . , µk )k∗,∞ = max

r∈[R]

Since each of our means has length at most

√

k
1 X ⊗r
µj
k
j=1

∗

.

(9)

d, we have that

k

kψ(µ1 , . . . , µk )k∗,∞ ≤ max

r∈[R]

1 X ⊗r
µj
k
j=1

13

∗

⊗r
R/2
≤ max max
.
√ kµ k∗ ≤ d
r∈[R] kµk≤ d

(10)

Using Claim 3.5, if |F| > N/η where


R/2

N = 1 + 2d

−R/4

/(d + 1)

D




D
3
d
3R/4
(2eR/d) · R log(5d) ,
≤ 1 + 2(d + 1)
≤ exp
4

we have that for at least 1 − η fraction of the Gaussian mixtures in F, there is another mixture
from F which is (d + 1)−R/4 close, as required.
Next we show that the closeness in moments implies closeness in the L2 distance. This follows
from fairly standard Fourier analytic techniques. We will first show that if the mean moments
are close, then the low-order Fourier coefficients are close. This will then imply that the Fourier
spectrum of the corresponding Gaussian mixtures f and f˜ are close.
Lemma 3.7. Suppose f (x), f˜(x) are the p.d.f. of G, G̃ which are both standard mixtures of k Gaussians√in d dimensions with means {µj : j ∈ [k]} and {µ̃j : j ∈ [k]} respectively that are both
ρ = d bounded. There exist universal constants c1 , c0 ≥ 1, such that for every ε ≤ exp(−c1 d) if
the following holds for R = c0 log(1/ε):
k
k
1 X ⊗r X
(µ̃j )⊗r
µj −
k

∀1 ≤ r ≤ R,

j=1

j=1

∗


≤ εr := ε

then kf − f˜k2 ≤ ε.

r
r
p
,
8πe log(1/ε)

(11)

Proof. We first show that if the moments are very close, the Fourier transform of the two distributions is very close. This translates to a bound on the L2 distance by Parseval’s identity.
Let g, h : Rd → R be defined by
h(x) =

k
k

X
1 X
δµ̃j (x) ,
δµj (x) −
k

2

g(x) = f (x) − f˜(x) = h(x) ∗ e−πkxk .

j=1

j=1

(12)

Since the Fourier transform of a convolution is the product of the Fourier transforms, we have
1 X −2πihζ,µj i  −πkζk2
∀ζ ∈ R , fb(ζ) =
e
·e
k
k

d

(13)

j=1

1 X −2πihµj ,ζi X −2πihµ̃j ,ζi 
b
h(ζ) =
,
e
e
−
k
k

k

j=1

j=1

2
gb(ζ) = b
h(ζ) · e−πkζk .

(14)

R
g (ζ)|2 dζ ≤ ε2 . We first note that the higher order Fourier coefficients
We will now show that Rd |b
of g do not contribute muchpto the Fourier mass of g because of the Gaussian tails. Let τ 2 =
4 log(1/ε). Since τ 2 ≥ (d + 2 d log(16/ε2 )) + 2 log(16/ε2 ))/(2π), using Lemma A.4
Z

kζk>τ

2

|b
g(ζ)| dζ =

Z

kζk>τ

2
|b
h(ζ)|2 e−2πkζk dζ ≤

14

Z

kζk>τ

2

4e−2πkζk dζ ≤

ε2
.
4

(15)

Now we upper bound |b
h(ζ)| for kζk < τ .

k
k ∞
1 X −2πihµj ,ζi
1 X X (−2πi)r
−2πihµ̃j ,ζi
b
|h(ζ)| =
(e
−e
) =
(hµj , ζir − hµ̃j , ζir )
k
k
r!
j=1 r=1

j=1

≤

∞
X
r=1

(2π)r kζkr 1
· ·
r!
k

k
X
j=1

µ⊗r
j −

k
X

(µ̃j )⊗r

j=1

∗

.

We claim that the injective norm above is at most kεr for all r ≥ 1. For r ≤ R, this follows
√
immediately from the assumption in (11). For r ≥ R, we use the fact that the means are ρ = d
bounded,
k
k
1 X ⊗r X
(µ̃j )⊗r
µj −
k
j=1

j=1

∗

≤2

max

µ∈{µj ,µ̃j }kj=1

≤ε·



kµkr2 ≤ 2dr/2 ≤ 2−r/2+1 (2d)r/2

r/2

r
R2
r
p
≤
ε
= εr ,
(8πe)2 log(1/ε)
8πe log(1/ε)

where the last line follows since R ≥ 16πe log(1/ε) and 2d ≤ log(1/ε). Hence, for kζk ≤ τ , we have
r

X (2πkζk)r
r
p
√
·ε
r!
2πr(r/e)r
8πe log(1/ε)
r
r=1

X 1  kζk r
X 1  2πekζk
r
ε
r
p
√
√
≤ε
≤ε
≤ ,
·
r
2
2πr
2πr 2τ
8πe log(1/ε)
r≥1
r≥1

|b
h(ζ)| ≤

∞
X
(2π)r kζkr

· εr ≤

since kζk ≤ τ . Finally, using this bound along with (15) we have
Z
Z
2
ε2
ε2 ε2
ε2
2
e−2πkζk dζ +
≤
+
≤ ε2 .
|b
g (ζ)| dζ ≤
4
4
4
4
kζk≤τ
Rd
Hence, by Parseval’s identity, the lemma follows.

The following lemma shows how to go from L2 distance to L1 distance √
using the CauchySchwartz lemma. Here we use the fact
p that all the means have length at most d. Hence, we can
focus on a ball of radius at most O( log(1/ε)), since both f, f˜ have negligible mass outside this
ball.
˜ two standard mixtures of Gaussians
Lemma 3.8. In the notation
√ above, suppose the p.d.f.s f, f of √
in d dimensions that are d-bounded (means having length ≤ d) satisfy kf − f˜k2 ≤ ε, for some
√
ε ≤ exp(−6d). Then, kf − f˜k1 ≤ 2 ε.
√
Proof. The means are all length at most d, and ε < 2−d . Let us define as before
τ2 =

√ p
p
p

1
d+2 d log(2/ε)+2 log(2/ε) , and γ = (maxkµk)+τ ≤ 2( d+ log(2/ε)) ≤ 4 log(2/ε).
2π
15

Let g = f − f˜ and let S = {x : kxk ≤ γ}. Using Gaussian concentration in Lemma A.4, we see
that both f, f˜ have negligible mass of at most ε/2 each outside S. Hence, using Cauchy-Schwartz
inequality and kf − f˜k2 ≤ ε,
sZ
sZ
Z
Z
p
g(x)2 dx ·
dx + ε ≤ ε Vol(S) + ε.
|g(x)|dx ≤ |g(x)|dx + ε ≤
Rd

S

S

S

Since S is a Euclidean ball of radius γ, by Stirling’s approximation (Fact A.6), the volume is
π d/2 γ d
≤
Vol(S) =
(d/2)!



2πeγ 2
d

d/2

≤



32πe log(2/ε)
d

d/2

≤ 2log(1/ε) =

1
,
ε

where the third inequality follows by raising both sides to the power 2/d and using the fact that
for α ≥ 6, 32πe(α + 1) ≤ 22α . This concludes the proof.
Proof of Theorem 3.2. The proof follows by a straightforward combination of Lemmas 3.6, 3.7,
and 3.8. As stated earlier, we will choose R = c0 log(1/ε), and constants c0 = 8πe, c1 = 36. First,
by Lemma 3.6 we have that for at least (1 − η) fraction of the Gaussian mixtures in F, there is
another mixture in the collection whose first R = c0 log(1/ε) mean moments differ by at most d−R/4
in symmetric injective tensor norm. To use Lemma 3.7 with ε′ = ε2 /4, we see that

r
r
ε2
p
min
min εr =
4 r∈[R] 8πe log(4/ε2 )
r∈[R]
√
ε2 −8πe√2 log(2/ε)
≥
·2
≥ 2−2 log(2/ε)−8πe 2 log(2/ε) ≥ 2−2πe log(1/ε) ≥ 2−R/4 ≥ (d + 1)−R/4 ,
4
as required, where for the first inequality we used that 2−α < 1/α for α > 0, and the second
inequality uses log(1/ε) ≥ c1 = 36. We complete the proof by applying Lemma 3.7 (with ε in
Lemma 3.7 taking the value ε2 /4) and Lemma 3.8.

4

√
√
Iterative Algorithms for min{Ω( log k), d} Separation

We now give a new iterative algorithm that estimates the means of a mixture of k spherical Gaussians up to √
arbitrary accuracy
δ > 0 in poly(d, k, log(1/δ)) time when the means have separation
√
of order Ω( log k) or Ω( d), when given coarse initializers. In all the bounds that follow, the
most interesting setting of parameters is when 1/δ is arbitrarily small compared to k, d (e.g.,
1/wmin ≤ poly(k) and δ = k−ω(1) , or when d = O(1) and δ = o(1)). For sake of exposition, we
will restrict our attention to the case when the standard deviations σi , and weights wi are known
for all i ∈ [k]. We believe that similar techniques can be used to handle unknown σi , wi as well
(see Remark 4.6). We also note that the results of this section are interesting even in the case of
uniform mixtures, i.e., wi = 1/k and σi = 1 for all i ∈ [k].
We will assume that we are given initializers that are inverse polynomially close in parameter
(0) (0)
(0)
distance. These initializers µ1 , µ2 , . . . , µk will be used to set up an “approximate” system of
non-linear equations that has sufficient “diagonal dominance” properties, and then use the Newton
method with the same initializers to solve it. In what follows, ρw and ρσ denote the aspect ratio
for the weights and variances respectively as defined in Section 2.
16

Theorem 4.1. There exist universal constants c, c0 > 0 such that the following holds. Suppose we
are given samples from a mixture of k spherical Gaussians G having parameters {(wj , µj , σj ) : j ∈
[k]}, where the weights and variances are known, satisfying
√
p
p
∀i 6= j ∈ [k], kµi − µj k2 ≥ c(σi + σj ) · min{ d + log(ρw ρσ ), log(ρσ /wmin )}
(16)
(0)

(0)

(0)

and suppose we are given initializers µ1 , µ2 , . . . , µk satisfying
∀j ∈ [k],

c0
1 (0)
kµ − µj k2 ≤
.
σj j
min{d, k}5/2

(17)

Then for any δ > 0, there is an iterative algorithm that runs in poly(ρ, d, 1/wmin , 1/δ) time (and
samples), and after T = O(log log(d/δ)) iterations recovers {µj : j ∈ [k]} up to δ relative error
(T )
(T )
w.h.p. i.e., finds {µj : j ∈ [k]} such that ∀j ∈ [k], we have kµj − µj k2 /σj ≤ δ.
√
√
For standard mixtures, (16) corresponds to a separation of order min{ log k, d}.
Firstly, we will assume without loss of generality that d ≤ k, since otherwise we can use a
PCA-based dimension-reduction result due to Vempala and Wang [45].
Theorem 4.2. Let {(wi , µi , σi ) : i ∈ [k]} be a mixture of k spherical Gaussians that is ρ-bounded,
and let wmin be the smallest mixing weight. Let µ′1 , µ′2 , . . . , µ′k be the projections onto the subspace
spanned by the top k singular vectors of sample matrix X ∈ Rd×N . For any ε > 0, with N ≤
−1
poly(d, ρ, wmin
, ε−1 ) samples we have with high probability
∀i ∈ [k], kµi − µ′i k2 ≤ ε.
The above theorem is essentially Corollary 3 in [45]. In [45], however, they take the subspace
spanned by the top max{k, log d} singular vectors, most likely due to an artifact of their analysis.
We give a different self-contained proof in Appendix C. We will abuse notation, and use {µi : i ∈ [k]}
to refer to the means in the dimension reduced space. Note that after dimension reduction, the
means are still well-separated for c′ > (c − 1) i.e.,

4.1

p
√
p
∀i, j ∈ [k], kµi − µj k ≥ c′ (σi + σj ) min{ d + log(ρw ρσ ), log(ρσ /wmin )}.

(18)

Description of the Non-linear Equations and Iterative Algorithm

For each component j ∈ [k] in G, we first define a region Sj around zj as follows. We will show in
Lemmas 4.8 and 4.9 that the total probability mass in Sj from other components is smaller than
the probability mass from the component j.
Definition 4.3 (Region Sj ). Given initializers z1 , z2 , . . . , zk ∈ Rd , define ebjℓ as the unit vector
along zℓ − zj , and let
√ p
p
Sj = {x ∈ Rd : ∀ℓ ∈ [k] |hx − zj , ebjℓ i| ≤ 4 log(ρσ /wmin )σj , and kx−zj k2 ≤ 4( d+ log(ρσ ρw ))σj }.
(19)
Based on those regions, we define the function F : Rkd → Rkd by

Z
k
1 X
(y − zj )gσi xi ,σi (y) dy,
wi
Fj (x) :=
wj σj
y∈Sj
i=1

17

(20)

for j = 1, . . . , k, where x = (x1 , . . . , xk ) ∈ Rkd . We also define x∗ = (µi /σi )ki=1 as the intended
solution. (Notice that for convenience, our variables correspond to “normalized means” µ/σ instead
of with the means µ directly.) The system of non-linear equations is then simply
F (x) = b

(21)

where b = F (x∗ ). This is a system of kd equations in kd unknowns. Obviously x∗ is a solution of
the system, and if we could find it we would be able to recover all the means, as desired.
Our algorithm basically just applies Newton’s method to solve (21) with initializers given by
(0)
(0)
xi = µi /σi for i ∈ [k]. To recall, Newton’s method uses the iterative update
x(t+1) = x(t) − (F ′ (x(t) ))−1 (b − F (x(t) )),
where F ′ (x(t) ) is the first derivative matrix (Jacobian) of F evaluated at x(t) . One issue, however, is
that we are not given b = F (x∗ ). Nevertheless, as we will show in Lemma 4.4, we can easily estimate
it to within any desired accuracy using a Monte Carlo approach based on the given samples (from
the Gaussian mixture corresponding to x∗ ). A related issue is that we do not have a closed-form
expression for F and F ′ , but again, we can easily approximate their evaluation at any point x and
to within any desired accuracy using a Monte Carlo approach (by generating samples from the
Gaussian mixture corresponding to x). The algorithm is given below in detail.

Iterative Algorithm for Amplifying Accuracy of Parameter Estimation
Input: Estimation accuracy δ > 0, N samples y (1) , . . . , y (N ) ∈ Rd from a mixture of well-separated
Gaussians G with parameters {(wj , µj , σj ) : j ∈ [k]}, the weights wj and variances σj2 , as well as
(0)

(0)

initializers µi for each i ∈ [k] such that kµi − µi k∞ ≤ ε0 σi .
Parameters: Set T = C log log(d/δ), for some sufficiently large constant
c0 d−5/2
√ C > 0, ε0 =
′
′
where c0 > 0 is a sufficiently small constant, and η1 , η2 , η3 = δwmin /(c dρσ ) where c > 0 is a
sufficiently small constant.
(T )
(T )
Output: Estimates (µi : i ∈ [k]) for each component i ∈ [k] such that kµi − µi k∞ ≤ δσi .
(T )

1. If δ ≥ ε0 , then we just output µi
(0)

2. Set xi

=

1 (0)
σi µi

(0)

= µi

for each i ∈ [k].

for each i ∈ [k].

3. Obtain using Lemma 4.4 an estimate b̃ of b up to accuracy η1 (in ℓ∞ norm). In more detail,
define the empirical average
∀j ∈ [k], b̃j =



X
1
I[y (ℓ) ∈ Sj ] y (ℓ) − zj .
wj σj N
ℓ∈[N ]

(This is the only place the given samples are used)
4. For t = 1 to T = O(log log(dk/δ)) steps do the following:
(a) Obtain using Lemma 4.4 an estimate F̃ (x(t) ) of F (x(t) ) at x(t) up to accuracy η2 (in ℓ∞
norm).
18

f′ (x(t) ) of F ′ (x(t) ) = ∇x F (x(t) ) up to accuracy η3
(b) Obtain using Lemma 4.5 an estimate F
(in ∞ → ∞ operator norm).

−1 

f′ (x(t) )
(c) Update x(t+1) = x(t) − F
b̃ − F̃ (x(t) ) .
(T )

5. Output µi

(T )

= σi xi

for each i ∈ [k].

The proof of the two approximation lemmas below is based on a rather standard Chernoff
argument; see Section 4.5.
Lemma 4.4 (Estimating Fj (x)). Suppose samples y (1) , y (2) , . . . , y (N ) are generated from a mixture
of k spherical Gaussians with parameters {(wi , σi xi , σi )}i∈[k] in d dimensions, and max{kxi k∞ , 1} ≤
ρ/σi ∀i ∈ [k]. There exists a constant C > 0 such that for any η > 0 the following holds for
2
N ≥ Cρ3 log( dk
γ )/(η wmin ) samples: with probability at least 1 − γ, for each j ∈ [k] the empirical
estimate for Fj (x) has error at most η i.e.,
Fj (x) −



X
1
I[y (ℓ) ∈ Sj ] y (ℓ) − zj
wj σj N
ℓ∈[N ]

∞

< η,

(22)

where Sj is defined in Definition 4.3.
We now state a similar lemma for the Jacobian F ′ : Rk·d → Rk·d , which by an easy calculation
is given for all i, j ∈ [k] by
Z
2πwi
(y − zj )(y − σi xi )T gσi xi ,σi (y) dy
∇xi Fj (x) =
wj σj σi y∈Sj
i.e., for all r, r0 ∈ [d],
∂Fj,r0 (x)
2πwi
=
∂xi (r)
wj σj σi

Z

y∈Sj

(y(r0 ) − zj (r0 ))(y(r1 ) − σi xi (r1 ))gσi xi ,σi (y) dy

(23)

Lemma 4.5 (Estimating ∇xi Fj (x)). Suppose we are given the parameters of a mixture of k spherical Gaussians {(wi , σi xi , σi )}i∈[k] in d dimensions, with max{kxi k∞ , 1} ≤ ρ/σi for each i ∈ [k],
and region Sj is defined as in Definition 4.3. There exists a constant C > 0 such that for any
2
(1) (2)
(N )
η > 0 the following holds for N ≥ Cρ4 d2 k2 log( dk
γ )/(η wmin ). Given N samples y , y , . . . , y
2
generated from spherical Gaussian with mean σi xi and variance σi /(2π), we have with probability
at least 1 − γ, for each j ∈ [k] the following empirical estimate for ∇xi Fj has error at most η i.e.,
^
∇xi Fj (x)−∇
xi Fj (x)

∞→∞

<

T


η
2πwi X
(ℓ)
(ℓ)
(ℓ)
^
y
−
σ
x
.
I[y
∈
S
]
y
−
z
, where ∇
F
(x)
:=
i i
j
j
xi j
k
wj σi σj N
ℓ∈[N ]

(24)
′
f
^
Furthermore, we have that the estimated second derivative F (x) = (∇xi Fj (x) : i, j ∈ [k]) satisfies
f′ (x)
F ′ (x) − F
≤ η.
∞→∞
19

Remark 4.6. Although for simplicity we focus here on the case that only the means are unknown,
we believe that our approach can be extended to the case that the weights and variances are also
unknown, and only coarse estimates for them are given. In order to handle that case, we need
to collect more statistics about the given samples, in addition to just the mean in each region.
Namely, using the samples restricted to Sj , we estimate the total probability mass in each Sj , i.e.,
(w)
(σ)
bj = Ey←G I[y ∈ Sj ], and the average squared Euclidean norm bj = d1 Ey←G I[y ∈ Sj ]kyk22 . We
now modify (21) by adding new unknowns (for the weights and variances), as well as new equations
(σ)
(w)
for bj and bj . This corresponds to k(d + 2) non-linear equations with k(d + 2) unknowns.

4.2

Convergence Analysis using the Newton method
(T )

We will now analyze the convergence of the Newton algorithm. We want each parameter xi ∈ Rd
to be close to x∗i in an appropriate norm (e.g., ℓ2 or ℓ∞ ). Hence, we will measure the convergence
and error of x = (xi : i ∈ [k]) to be measured in ℓ∞ norm.
Definition 4.7 (Neighborhood). Consider a mixture of Gaussians with parameters ((µi , σi , wi ) :
i ∈ [k]), and let (xi : i ∈ [k]) ∈ Rkd be the corresponding parameters of the non-linear system
F (x) = b. The neighborhood set
N = {(xi : i ∈ [k]) ∈ Rkd | ∀i ∈ [k], kxi − x∗i k∞ < ε0 = c0 d−5/2 },
is the set of values of the variables that are close to the true values of the variables given by
x∗i = µσii ∀i ∈ [k], and c0 > 0 is an appropriately large universal constant given in Theorem 4.1.
We will now show the convergence of the Newton method by showing diagonal dominance
properties of the non-linear system given in Lemma 2.5. This diagonal dominance arises from the
separation between the means of the components. Lemmas 4.8 and 4.9 show that most of the
probability mass from jth component around µj is confined to Sj , while the other components of
G are far enough from zj , that they do not contribute much ℓ1 mass in total to Sj .
The following lemma lower bounds the contribution to region Sj from the jth component.
Lemma 4.8. In the notation of Theorem 4.1, for all j ∈ [k], we have
Z
1
.
gµj ,σj (y) dy ≥ 1 −
8πd
Sj
Z
1
1
∀r1 ∈ [d],
(y(r1 ) − µj (r1 )) gµj ,σj (y) dy ≤
.
σj Sj
8πd
Z
σj2
1
1
I[r1 = r2 ] ≤
.
(y(r1 ) − µj (r1 )) (y(r2 ) − µj (r2 )) gµj ,σj (y) dy −
∀r1 , r2 ∈ [d], 2
2π
8πd
σj Sj

(25)
(26)
(27)

The proof of the above lemma follows from concentration bounds for multivariate Gaussians.
The following lemma upper bounds the contribution in Sj from the other components.
Lemma 4.9 (Leakage from other components). In the notation of Theorem 4.1, and for a component j ∈ [k] in G, the contribution in Sj from the other components is small, namely, for all

20

j ∈ [k]
X

wi

i∈[k],i6=j

∀r1 ∈ [d],
∀r1 , r2 ∈ [d],

X

i∈[k],i6=j

X

i∈[k],i6=j

Z

Sj

gµi ,σi (y) dy <

wj
.
16πd

Z
wj
wi kµi − µj k2
·
.
|y(r1 ) − µi (r1 )| gµi ,σi (y) dy <
σi
σi
16πdρσ
Sj
Z
wj
wi
|y(r1 ) − µi (r1 )| |y(r2 ) − µi (r2 )| gµi ,σi (y) dy <
.
σi σj Sj
16πdρσ

(28)
(29)
(30)

The above lemma is the more technical of the two, and it is crucial in showing diagonal√dominance
√of the Jacobian. The proof of the above lemma is very different for separation of order log k
and d – hence we will handle this separately in Sections 4.4.1 and 4.4.2 respectively.
We now show the convergence of Newton algorithm assuming the above two lemmas. Theorem 4.1 follows in a straightforward manner from the guarantees of the Newton algorithm. We
mainly need to show that kF ′ kkF ′′ kε0 < 1/2.
We now prove that the function (F ′ (x))−1 has bounded operator norm using the diagonal
dominance properties of F .
Lemma 4.10. For any point x ∈ N , the operator F ′ : Rd·k → Rd·k satisfies k(F ′ (x))−1 k∞→∞ ≤ 4.
Proof. We will divide the matrix F ′ into k × k = k2 blocks of size d × d each, and show that the
matrix satisfies the required diagonal dominance property. Let us first consider the diagonal blocks
i.e., we have from (23) for each j ∈ [k]
Z
2π
(y − zj )(y − σj xj )T gσj xj ,σj (y) dy
∇xj Fj (x) = 2
σj y∈Sj
Z
Z
2π(σj xj − zj )
2π
(y − σj xj )(y − σj xj )T gσj xj ,σj (y) dy +
(y − σj xj )T gσj xj ,σj (y) dy
= 2
σj y∈Sj
σj2
y∈Sj
Consider
a mixture of Gaussians where the jth component has mean σj xj and standard deviation
√
σj / 2π.
It satisfies the required separation conditions since kσi xi − µi k2 ≤ ε0 σi . Applying
Lemma 4.8,
Z
2π
(y − σj xj )(y − σj xj )T gσj xj ,σj (y) dy
∇xj Fj (x) = Id×d − 2
σj y∈S
/ j
Z
2π(σj xj − zj )
(y − σj xj )T
−
gσj xj ,σj (y) dy
σj
σj
y ∈S
/ j
= Id×d − E,

where E ∈ Rd×d satisfies from (27) and (26) that
1
1
1
∀r1 , r2 ∈ [d], |E(r1 , r2 )| ≤
+ ε0 ·
≤
4d
4d
2d
X
1
∀r1 ∈ [d],
|E(r1 , r2 )| ≤ .
2
r2 ∈[d]

21

(31)

Using a similar calculation, we see that for the off-diagonal blocks Mji = ∇xi Fj (x),
Z
2πwi
(y − zj )(y − σi xi )T gσi xi ,σi (y) dy
Mji = ∇xi Fj (x) =
wj σi σj y∈Sj
Z
2πwi
(y − σi xi )(y − σi xi )T
=
gσi xi ,σi (y) dy
wj y∈Sj
σi σj
Z
(y − σi xi )T
2πwi (σi xi − zj )
+
gσi xi ,σi (y) dy
wj σj
σi
y∈Sj
Also, kσi xi − zj k∞ ≤ kµi − µj k∞ + σi + σj ≤ 2kµi − µj k2 . For each r1 , r2 ∈ [d]
Z
2πwi ρσ
|Mji (r1 , r2 )| ≤
|y(r1 ) − σi xi (r1 )| |y(r2 ) − σi xi (r2 )| gσi xi ,σi (y) dy
wj σi2 y∈Sj
Z
|y(r2 ) − σi xi (r2 )|
4πwi ρσ kµi − µj k2
gσi xi ,σi (y) dy.
+
wj
σi
σi
y∈Sj
Summing over all i 6= j, and using the bounds in Lemma 4.9,
X

i∈[k],i6=j

|Mji (r1 , r2 )| ≤

1
4d

d X
X
1
|Mji (r1 , r2 )| ≤
4

(32)

r2 =1 i6=j

Hence using Lemma 2.5 for diagonally dominant matrices, we get from (31) and (32)
k(F ′ (x))−1 k∞→∞ ≤

1
1−

1
2

−

1
4

≤ 4.

We now prove that the function F ′ (x) is locally Lipschitz by upper bounding the second derivative operator.
Lemma 4.11. There exists a universal constant c′ > 0 such that the derivative F ′ is locally LLipschitz in the neighborhood N for L ≤ c′ d5/2 .
Proof. We proceed by showing that the operator norm of F ′′ : Rdk×dk → Rdk is upper bounded by
L at any point x ∈ N . We first write down expressions for F ′′ , and then prove that the operator
norm of F ′′ is bounded. We first observe that for each j, i1 , i2 ∈ [k],
∀r0 , r1 , r2 ∈ [d],

∂ 2 Fj,r0 (x)
= 0 if i1 6= i2 .
∂xi1 (r1 )∂xi2 (r2 )

Hence the second derivatives are non-zero only for diagonal blocks i1 = i2 . For each j, i ∈ [k] the

22

second derivatives for ∀r0 , r1 , r2 ∈ [d] are given by
Z

 (y(r ) − σ x (r ))(y(r ) − σ x (r ))

1
i i 1
2
i i 2
(y(r0 ) − zj (r0 ))
−
I[r
=
r
]
gσi xi ,σi (y) dy
1
2
σi2
y∈Sj
Z

 (y(r ) − σ x (r ))(y(r ) − σ x (r ))
∂ 2 Fj,r0 (x)
4πwi
ky − zj k2
1
i i 1
2
i i 2
+
1
gσi xi ,σi (y) dy.
≤
max
·
∂xi (r1 )∂xi (r2 )
wj y∈Sj
σj
σi2
y∈Sj
√ Z
 (y(r ) − σ x (r ))(y(r ) − σ x (r ))

8πwi d
1
i i 1
2
i i 2
≤
+
1
gσi xi ,σi (y) dy.
wj
σi2
y∈Sj
∂ 2 Fj,r0 (x)
4πwi
=
∂xi (r1 )∂xi (r2 )
wj σj

A simple bound on the operator norm of F ′′ : Rkd × Rkd → Rkd is given by summing over all i ∈ [k]
and r1 , r2 ∈ [d]. Consider
√ a mixture of Gaussians where the jth component has mean σj xj and
standard deviation σj / 2π. It satisfies the required separation conditions since kσi xi −µi k2 ≤ ε0 σi .
Applying Lemma 4.9 we get
kF ′′ (x)k∞,∞→∞ ≤
X X

i6=j r1 ,r2 ∈[d]

X

r1 ,r2 ∈[d]

X

X

∂ 2 Fj,r0 (x)
∂xi (r1 )∂xi (r2 )

i∈[k] r1 ,r2 ∈[d]

∂ 2 Fj,r0 (x)
wj
16πd5/2
≤
≤ d3/2 , and
·
∂xi (r1 )∂xi (r2 )
wj
16πd
∂ 2 Fj,r0 (x)
≤ 8πd5/2
∂xj (r1 )∂xj (r2 )

Z

y∈Rd

+ 8πd5/2

Z

(y(r1 ) − σj xj (r1 ))(y(r2 ) − σj xj (r2 ))
gσj xj ,σj (y)
σj2

y∈Rd

gσj xj ,σj (y) dy


1
+
1
≤ 15πd5/2
π2
≤ 16πd5/2 .
≤ 8πd5/2

∀x ∈ N , kF ′′ (x)k∞,∞→∞

Hence, applying Lemma B.1 with F ′ in the open set K = N , the lemma follows.
(0)

Proof of Theorem 4.1. We first note that we have from (17), for each i ∈ [k] that kx∗i − xi k∞ ≤
√
(0)
(0)
/σi ≤ ε0 . We will use the algorithm to find kx(T ) − x∗ k∞ < δ/ d (use
kx∗i − xi k2 ≤ kµi − µi k2√
the algorithm with δ′ := δ/ d). The proof follows in a straightforward manner from Corollary B.4.
Lemma 4.11 shows that F ′ is locally L-Lipschitz for L = c′ d5/2 . Together with Lemma 4.10 we have
−3
)
ε0 Lk(F ′ )−1 k∞→∞ ≤ 1/2 for our choice of ε0 . Also from Lemma 4.4, by using N = O(δ−2 ρ6 k6 wmin
. It is easy to check that B = maxy kF ′ (y)k ≤
samples (and d ≤ k), η1 + η2 ≤ 4√dk(F ′ )δ−1 k
maxi σi
wmin mini σj

∞→∞

≤ 4ρσ /wmin . Hence, similarly from Lemma 4.5 η3 ≤

δ
√
4 dBk(F ′ )−1 k2∞→∞

as required.

Hence, from Corollary B.4, T = O(log log(d/δ)) iterations of the Newton’s method suffices. Further,
each iteration mainly involves inverting a dk × dk matrix which is polynomial time in d, k.

23

4.3

Bounding Leakage from Own Component: Proof of Lemma 4.8

Lemma 4.12. Let vb0 ∈ Rd be a unit vector. Then for any coordinates r1 , r2 ∈ [d] we have for any
t≥3
Z
1
1
exp(−πkyk22 /σ 2 ) dy < 16π
exp(−t2 )
(33)
σ d y:hy,bv0 i≥tσ
Z
1
1
exp(−t2 )σ
(34)
|y(r1 )| exp(−πkyk22 /σ 2 ) dy < 16π
σ d y:hy,bv0 i≥tσ
Z
1
1
exp(−t2 )σ 2
(35)
|y(r1 )y(r2 )| exp(−πkyk22 /σ 2 ) dy < 16π
σ d y:hy,bv0 i≥tσ
Proof. The first equation (33) follows easily from the rotational invariance of a spherical Gaussian.
Suppose g0 = hy, vb0 i, then g0 ∼ N (0, σ 2 /(2π)). Hence
Z
Z
√ 
1
1
2
2
2
2
2πt < exp(−t2 ).
exp(−πkyk
/σ
)
dy
=
exp(−πg
/σ
)
dg
≤
Φ
0
0,1
2
0
σ d y:hy,bv0 i≥tσ
σ g0 ≥tσ
We will now prove (35); the proof of (34) follows the same argument. We first observe that
using the rotational invariance, it suffices to consider the span of vb0 , ebr1 , ebr2 .
Consider the case r1 6= r2 . Let vb1 , vb2 be the unit vectors along ebr1 − hb
er1 , vb0 ib
v0 and ebr2 −
hb
er2 , vb0 ib
v0 −hb
er2 , vb1 ib
v1 respectively. If gℓ = hy, vbℓ i for ℓ = 0, 1, 2, then {gℓ } are mutually independent
and distributed as N (0, σ 2 /(2π)). Let α1 = hb
v0 , ebr1 i, α2 = hb
v0 , ebr2 i and β2 = hb
v1 , ebr2 i. Then,
q
q
y(r1 ) = α1 g0 + 1 − α21 · g1 , and y(r2 ) = α2 g0 + β2 g1 + 1 − α21 − β22 · g2 ,
|y(r1 )| ≤ |g0 | + |g1 |, and |y(r1 )y(r2 )| ≤ (|g0 | + |g1 |)(|g0 | + |g1 | + |g2 |) ≤ 3(|g0 |2 + |g1 |2 + |g2 |2 ).

Hence,
Z
Z

3
1
2
|y(r1 )y(r2 )| exp(−πkyk2 ) dy ≤ d
|g0 |2 + |g1 |2 + |g2 |2 exp(−πkyk22 ) dy
d
σ y:hy,bv0 i≥tσ
σ y:hy,bv0 i≥tσ
Z
Z Z
3
≤ 3
(|g0 |2 + |g1 |2 + |g2 |2 ) · exp(−πg22 ) dg2 · exp(−πg12 ) dg1 · exp(−πg02 ) dg0
σ g0 >tσ g1 g2
Z
Z
Z
3
6
1
2
2
2
2
≤
|g0 | exp(−πg0 ) dg0 +
exp(−πg02 ) dg0 (by symmetry)
|g1 | exp(−πg1 ) dg1 ·
σ g0 >tσ
σ g1
σ g0 >tσ
≤ 3 exp(−2t2 )σ 2 + 6σ 2 exp(−2t2 ) ≤ 9σ 2 exp(−2t2 ),

where the last line follows from truncated moments of a normal variable with variance σ 2 /(2π)
(Lemma A.2). Using the fact that exp(t2 ) > 9 · 16π for t ≥ 3, the lemma follows. The proof for
r1 = r2 , and for (34) follows from an identical argument.
Proof of Lemma 4.8. We will consider the loss from the region outside of Sj . The region Sj is
defined by k linear inequalities, one for each of the components ℓ ∈ [k], and a constraint that points
lie in an ℓ2 ball around zj . Consider y drawn from gµj ,σj . Since kzj − µj k2 ≤ σj and Lemma A.4
(applied with t = 6d),
√
√
P[ky − zj k2 ≥ 4 dσj ] ≤ P[ky − µj k2 ≥ 3 dσj ] ≤ exp(−6d) ≤
24

1
.
16πd

When y is drawn from gµj ,σj , from Lemma 4.12, each of the k linear constraints are satisfied with
high probability i.e., for a fixed ℓ ∈ [k],
h
i
h
i w2
p
p
1 wmin
,
P |hy − zj , ebjℓ i| > 4σj log(1/wmin ) ≤ P |hy − µj , ebjℓ i| > 3σj log(1/wmin ) ≤ min ≤ ·
16π
k 16π
(36)
where the last step follows since wmin ≤ 1/k. Hence, performing a union bound over the k components and the ℓ2 constraint, and using d ≤ k, we get
Z
1
1
1
gµj ,σj (y) dy ≥ 1 −
−
≥1−
.
16πd
16πd
8πd
Sj
Z

To prove (26), ∀r1 ∈ [d] we have
Z
Z
(y(r1 ) − µj (r1 )) gµj ,σj (y) dy −
(y(r1 ) − µj (r1 )) gµj ,σj (y) dy
(y(r1 ) − µj (r1 )) gµj ,σj (y) dy =
y ∈S
/ j
Sj
Rd
Z
=0−
(y(r1 ) − µj (r1 )) gµj ,σj (y) dy
y ∈S
/ j

p
Any point y ∈
/√Sj implies that one of the k linear constraints hy − µj , ebji i > 3 log(1/wmin )σj , or
ky − µj k2 > 3 dσj is true. Again, for the ℓ2 ball constraint, from Lemma A.5 with q ≤ 2,
Z

(y(r1 ) − µj (r1 ))
· gµj ,σj (y) dy ≤
σj

√
ky−µj k2 >3 dσj

Z

ky − µj k2
1
· gµj ,σj (y) dy ≤
.
σj
16πd

√
ky−µj k2 >3 dσj

p
Let Zji = {y : hy − µj , ebji i > 3σj log(ρσ /wmin )} be the set y that do not satisfy the constraint
along ebji . Since kzj − µj k2 ≤ σj , we have that the total loss from those y ∈ Zji is
Z

y∈Zji

Z

(y(r1 ) − µj (r1 )) gµj ,σj (y) dy ≤

y∈Zji

≤ σj



|y(r1 ) − µj (r1 )| gµj ,σj (y) dy

2
wmin
16πρ2σ



,

p
by applying Lemma 4.12 with the Gaussian (y − µj ), and t = 3 log(ρσ /wmin ). Hence, we have
Z

Sj

(y(r1 ) − µj (r1 )) gµj ,σj (y) dy ≤

Z

y ∈S
/ j

≤ σj



(y(r1 ) − µj (r1 )) gµj ,σj (y) dy

2
wmin
16πρ2σ



k+

 1 
σj
,
≤ σj
16πd
8πd

since wmin ≤ 1/k. The proof of the last equation follows in an identical fashion. From (35) of
Lemma 4.12, we have
1
∀r1 , r2 ∈ [d], 2
σj

Z

y ∈S
/ j



1
y(r1 ) − µj (r1 ) y(r2 ) − µj (r2 ) gµj ,σj (y) dy ≤
.
8πd
25

Further,
Z

σj2
I[r1 = r2 ]
(y(r1 ) − µj (r1 )) (y(r2 ) − µj (r2 )) gµj ,σj (y) dy =
2π
Rd

∀r1 , r2 ∈ [d],
Hence,

4.4

1
σj2

Z

Sj

(y(r1 ) − µj (r1 )) (y(r2 ) − µj (r2 )) gµj ,σj (y) dy −

1
1
I[r1 = r2 ] ≤
.
2π
8πd

Bounding Leakage from Other Components: Proof of Lemma 4.9

√
Since
our algorithm works when we either have separation of order log kpor separation of√order
√
d, we have two different proofs for Lemma 4.9 depending on whether log(1/wmin ) ≤ d or
not3 .
4.4.1

Separation of Order

√
log k

In this case, we have a separation of
p
∀i 6= j ∈ [k], kµi − µj k2 ≥ c(σi + σj ) log(ρσ /wmin ).

Let ebji be the unit vector along zi − zj . We will now show that every point y ∈ Sj , is far from
zi along the direction ebji .
p
hzi − zj , ebji i = kzj − zi k2 ≥ kµi − µj k2 − σi − σj ≥ 5 log(ρσ /wmin )(σi + σj ) + 21 kµi − µj k2
p
Since y ∈ Sj , |hy − zj , ebji i| ≤ 4σj log(ρσ /wmin ), and
p
hzi − y, ebji i ≥ hzi − zj , ebji i − |hy − zj , ebji i| ≥ 5σi log(ρσ /wmin ) + 21 kµi − µj k2
p
hµi − y, ebji i ≥ hzi − y, ebji i − kµi − zi k2 ≥ 4σi log(ρσ /wmin ) + 12 kµi − µj k2

We will now use Lemma 4.12 for each component i ∈ [k] with mean zero Gaussian y − µi and
1
2
t2 = 16 log(ρσ /wmin ) + 8πσ
2 kµj − µi k2 . We first prove (28); using equation (33) with Gaussian
i

(y − µi ) of variance σi2 /(2π),
Z
 kµ − µ k2  w2
w2
i
j 2
< min .
gµi ,σi (y) dy ≤ wi exp(−t2 ) < wi · min exp −
wi
2
16π
16π
16σi
Sj
Z
X
wj
.
gµi ,σi (y) dy <
Hence,
wi
16πd
Sj
i∈[k],i6=j

For (29), we use (34) similarly with (y − µi ) and variance σi2 /(2π) and t as above:
Z
 kµ − µ k2 
|y(r1 ) − µi (r1 )|
w2
i
j 2
gµi ,σi (y) dy ≤ wi · min exp −
wi
.
2
σ
16πρ
16σ
i
σ
Sj
i
X wi kµi − µj k2 Z
wj
|y(r1 ) − µi (r1 )| gµi ,σi (y) dy <
·
.
σi
σi
16πdρ
σ
Sj
i∈[k],i6=j

3

More accurately whether

p
√
p
log(ρσ /wmin ) ≤ d + log(ρw ρσ ), where ρw ≥ 1/wmin .

26

The proof of (30) follows in an identical fashion from (35) for r1 , r2 ∈ [d]
Z
2
wj wi
wi σi wmin
wi
|y(r1 ) − µi (r1 )| |y(r2 ) − µi (r2 )| gµi ,σi (y) dy ≤
·
≤
2
σi σj Sj
σj 16πρσ
16πdρσ
Z
X
wj
wi
.
|y(r1 ) − µi (r1 )| |y(r2 ) − µi (r2 )| gµi ,σi (y) dy <
σi σj Sj
16πdρσ
i∈[k],i6=j

4.4.2

Separation of Order

√
d

√
The proof of Lemma 4.9 uses the following useful lemma that shows that for any point within σj d
distance of µj , the total probability mass from the other Gaussian components is negligible. Note
that in the following lemma, k can be arbitrary large in terms of d.
Lemma 4.13. There exists a universal constant C0 > 1 such that for any C ≥ C0 , if G is a mixture
of k spherical Gaussians {(wj , µj , σj ) : j ∈ [k]} in Rd satisfying
√

p
∀i 6= j ∈ [k], kµi − µj k2 ≥ C(σi + σj )
d + log(ρw ρσ ) ,
(37)

√
then for every component j ∈ [k] and ∀x∗ : kx∗ − µj k ≤ 4σj d and ∀0 ≤ m ≤ 2,
X
wi gµi ,σi (x∗ ) < exp(−C 2 d/8) · wj gj (x).
X

i∈[k],i6=j

wi



(38)

i∈[k],i6=j

kx − µi km
σim



· gµi ,σi (x∗ ) ≤

exp(−C 2 d/8)
· wj gj (x∗ ),
ρm
σ

(39)

where ρσ = maxi σi /(mini σi ) and ρw = maxi wi /(mini wi ).
p
√
We now prove Lemma 4.9 under separation kµi − µj k2 ≥ c(σi + σj )( d + log(ρσ ρw )) assuming
the above lemma.
Proof for Lemma 4.9. Our proof will follow by
Lemma 4.13, and integrating over all
p
√ applying
y ∈ Sj . For any point y ∈ Sj , ky − µj k2 ≤ 3σj ( d + log(ρw ρσ )), and the separation of the means
satisfies (37). To prove (25) we get from (38),
X
wi gµi ,σi (y) ≤ wj exp(−C 2 d/8) · gµj ,σj (y).
i∈[k],i6=j

Integrating over Sj ,

X

i∈[k],i6=j

wi

Z

Sj

gµi ,σi (y) dy ≤ wj exp(−C 2 d/8)

Z

Rd

gµj ,σj (y) dy ≤

wj
.
16πd

To prove (27), we get from (39) that for each y ∈ Sj and r1 , r2 ∈ [d],
X

i∈[k]\{j}

wi
|y(r1 ) − µi (r1 )| |y(r2 ) − µi (r2 )| gµi ,σi (y) ≤
σi σj
≤
27

X

i∈[k]\{j}

wi ρσ ky − µi k22
gµi ,σi (y)
wj
σi2

wj exp(−C 2 d/8)
· gµj ,σj (y).
ρσ

Integrating over all y ∈ Sj ,
Z
Z
X
wi
wj
wj
2
exp(−C d/8)·
gµj ,σj (y) ≤
|y(r1 ) − µi (r1 )| |y(r2 ) − µi (r2 )| gµi ,σi (y) dy ≤
σi σj Sj
ρσ
16πdρσ
Sj

i∈[k],i6=j

√
√
To prove (26), we first note that for each y ∈ Sj , kµi − µj k2 ≤ kµi − yk + ( d + log ρw )σj ≤
2kµi − yk. Hence, by applying (39) in Lemma 4.13, the following inequality holds:
X

i∈[k],i6=j

 ky − µ k 2
X
wi kµi − µj k2 ky − µi k2
i
gµi ,σi (y) ≤ 2
wi
·gµi ,σi (y) ≤ 2 exp(−C 2 d/8)ρ−2
σ ·wj gj (y).
σi σj
σi
i∈[k],i6=j

(40)

Hence (26) follows from a similar argument as before.
The proof of Lemma 4.13 proceeds by using a packing-style argument. Roughly speaking, in
the uniform case when all the variances are roughly equal, the separation condition can be used to
establish an upper bound on the number of other Gaussian means that are at a distance of r from
a certain mean. We now present the proof for the case when the variances are roughly equal, and
this will also be important for the general case.
Lemma 4.14. In the notation of Lemma 4.13, let us denote for convenience, the jth√component
by (w0 , µ0 , σ0 ). Then, the total probability mass at any point x∗ s.t. kx∗ − µ0 k ≤ 4 dσ0 , from
components (w1
, µ1 , σ1 ), . . . , (wi , µi , σi ), . . . , (wk′ , µk′ , σk′ ) satisfying ∀i ∈ [k′ ], σ0 ≤ σi ≤ 2σ0 , and
√
√
d + log ρw is given by
kµi − µ0 k2 ≥ C
X

wi gi (x∗ ) < w0 σ0−d exp(−C 2 d/4).

(41)

i∈[k ′ ]

Furthermore, we have

P

1/2 ∗
i∈[k ′ ] wi gi (x )

< w0 σ0−d exp(−C 2 d/4).

Proof. We can assume without loss of generality that µ0 = 0, σ0 = 1. Suppose ri = kµi − xk2 .
Considering 1 ≤ σi ≤ 2 and gi (x) = σi−d exp(−πkx − µi k22 /σi2 ), the lemma will follow, if we prove
the following inequality for any C ≥ C0 ≥ 16:
X

i∈[k ′ ]

Let Bi = {y : ky − µi k2 ≤

√


wi exp − π2 · ri2 < w0 exp(−C 2 d/4),

(42)

d}. From the separation conditions, we have that

• the balls B1 , . . . , Bk′ are pairwise disjoint, and

√
√
• the balls are far from the origin i.e, ∀y ∈ Bi , kyk ≥ C( d + log ρw ).

√
We first relate the p.d.f. value gi (x) to the Gaussian measure of a ball of radius d around µi .
The volume of the ball Vol(Bi ) ≥ 1 and for every y ∈ Bi , the√separation conditions and triangle
inequality imply that kyk ≥ kµi − x∗ k − ky − µi k − kx∗ k ≥ ri / 2. Hence,
 Z


πri2
exp −πkyk2 dy.
(43)
≤
exp −
2
y∈Bi
28

By using the disjointness of balls Bi , we get that


Z
XZ
X
X


πri2
2
exp −πkyk2 dy
exp −πkyk dy ≤ wmax
≤
wi
wi exp −
2
y∈Bi
i∈[k ′ ] y∈Bi
i∈[k ′ ]
i∈[k ′ ]
Z

≤ wmax
exp −πkyk2 dy
√ √
y:kyk≥C( d+ log ρw )

≤ ρw w0 ·

1
exp(−C 2 d/4) ≤ w0 exp(−C 2 d/4),
ρw

where the last line follows since a Gaussian random variable in d dimensions with mean
√ 0 and unit
variance in each direction has measure at most exp(−s2 /2) outside a ball of radius ( d + s) (see
Lemma A.4).
We now proceed to the proof of Lemma 4.13.
Proof of Lemma 4.13. We may again assume without
√ √ loss of generality that µj = 0 and σj = 1 (by
shifting the origin and scaling). Hence kx∗ k2 ≤ d/ π. We will divide the components i ∈ [k]\{j}
depending on their standard deviation σi into buckets I0 , I1 , . . . , Is where s ≤ ⌈log(maxi σi /σj )⌉ as
follows:
∀q ∈ [s] Iq = {i ∈ [k] : 2q−1 σj < σi ≤ 2q σj }.

I0 = {i ∈ [k] \ {j} : σi ≤ σj },

Let us first consider the components in the bucket I0 , and suppose we scale so√that σj = 1.
As before let ri = kx∗ − µi k2 . We first note that if σi ≤ σj , since kx∗ − µi k2 ≥ C d(σi + σj ), a
simple calculation shows that gµi ,σi (x∗ ) ≤ gµi ,σj (x∗ ). Hence, by applying Lemma 4.14 to I0 with a
uniform variance of σj2 = 1 for all Gaussians, we see that
X

i∈I0

wi gµi ,σi (x∗ ) ≤

X

i∈I0

wi gµi ,σj (x∗ ) ≤ wj exp(−C 2 d/4).

(44)

Consider any bucket Iq . We will scale the points down by 2q−1 σj so that y ′ = y/(2q−1 σj ), and let
∀i ∈ Iq ri′ = ri /(2q−1 σj ). Again from Lemma 4.14 we have that
X

i∈Iq

Hence,

X

i∈Iq

wi exp(−πri2 /σi2 ) ≤

wi gi (x) ≤

1
2(q−1)d σjd

X

(45)

wi exp(−πri2 /σi2 ) ≤ wj 2−(q−1)d exp(−C 2 d/4).

(46)

i∈Iq

X

i∈Iq

2

wi exp(−πri′ /4) ≤ wj exp(−C 2 d/4).

P
Hence, by summing up over all buckets (using (44) and (46)), the total contribution i∈[k]\{j} wi gi (x) ≤
4wj exp(−C d /4) ≤ exp(−C 2 d/8) · wj gj (x).
√
The final equation (39) follows from separation conditions, since kx − µi k > C dσi , hence
exp(πkx − µi k2 /2σi2 ) > kx − µi km /σim for some sufficiently large constant C = C(m) ≥ 1. Hence,
by using an identical argument with the furthermore part of Lemma 4.14, it follows.

29

4.5

Sampling Errors

Lemma 4.15 (Error estimates for Gaussians). Let S ⊂ {x ∈ Rd : kxk2 ≤ ρ} be any region,
and suppose samples x(1) , x(2) , . . . , x(N ) are generated from a mixture of k spherical Gaussians with
parameters {(wi , µi , σi )}i∈[k] in d dimensions, with ∀i, kµi k, σi ≤ ρ. There exists a constant C > 0
such that for any ε > 0, with N ≥ C ε−2 ρ2 log d log(1/γ) samples, we have for all r, r ′ ∈ [d] with
probability at least 1 − γ
Z

X

(x(r) − µi (r))gµi ,σi (x)dx −

x∈S i∈[k]

Z

X


1 X  (ℓ)
x (r) − µi (r) I[x(ℓ) ∈ S] < ε.
N

(47)

ℓ∈[N ]

(x(r) − µi (r))(x(r ′ ) − µi (r ′ ))T gµi ,σi (x)dx

x∈S i∈[k]

−


T
1 X  (ℓ)
x (r) − µi (r) x(ℓ) (r ′ ) − µi (r ′ ) I[x(ℓ) ∈ S] < ε.
N

(48)

ℓ∈[N ]

Proof. Fix an element r ∈ [d]. Each term ℓ ∈ [N ] in the sum corresponds to an i.i.d random
variable
1 P
(ℓ)
(ℓ)
Zℓ = x (r) − µi (r) I[x ∈ S]. We are interested in the deviation of the sum Z = N ℓ∈[N ] Zℓ .
R
P
Firstly, E[Z] = x∈S i∈[k] (x(r) − µi (r))gµi ,σi (x)dx. Further, each of the i.i.d r.v.s has value

|Zℓ − E Zℓ | ≤ |x(ℓ) (r) − µi (r)| + ρ. Hence, |Zℓ | > (2ρ + t maxi σi ) with probability O exp(−t2 /2) .
Hence, by using standard sub-gaussian tail inequalities, we get


ε2 N
P[|Z − E Z| > ε] < exp −
(2ρ + maxi σi )2

Hence, to union bound over all d events N = O ε−2 ρ2 log d log(1/γ) suffices.
A similar proof also works for the second equation (48).

P
Proof of Lemma 4.4. Let Z = wj σ1j N ℓ∈[N ] I[y (ℓ) ∈ Sj ] y (ℓ) − zj . We see that
Z
k
1 X
E[Z] =
Ij (y)(y − zj )gσi xi ,σi (y) dy = Fj (x).
wi
wj σj
y∈Rd
i=1

Further, we can write Z = Z1 + Z2 , where
Z1 =



X
1
I[y (ℓ) ∈ Sj ] y (ℓ) − σi xi ,
wj σj N
ℓ∈[N ]

X
1
(σi xi − zj )
I[y (ℓ) ∈ Sj ].
Z2 =
wj σj N
ℓ∈[N ]

By applying Lemma 4.15, we get
kZ − E[Z]k∞ ≤ kZ1 − E Z1 k∞ + kZ2 − E Z2 k∞ ≤

30

η η
+ = η.
2 2

T

P
Proof of Lemma 4.5. Let Z = wj σwi σi j N ℓ∈[N ] I[y (ℓ) ∈ Sj ] y (ℓ) − zj y (ℓ) − σi xi , where the samples are drawn just from the spherical Gaussian with mean σi xi and variance σi2 . Hence,
Z
wi
E[Z] =
Ij (y)(y − zj )(y − σi xi )T gσi xi ,σi (y) dy = Fj (x).
wj σj σi y∈Rd
Also Z = Z1 + Z2 where
T

X
wi
and
I[y (ℓ) ∈ Sj ] (σi xi − zj ) y (ℓ) − σi xi
wj σi σj N
ℓ∈[N ]
T


X
wi
I[y (ℓ) ∈ Sj ] y (ℓ) − σi xi y (ℓ) − σi xi .
Z2 =
wj σi σj N

Z1 =

ℓ∈[N ]

The kZ − E Zk∞→∞ is the sum of the absolute values of the d entries in a row. From Lemma 4.15,
η
η
kZ − E Zk∞→∞ = kZ1 − E Z1 k∞→∞ + kZ2 − E Z2 k∞→∞ ≤
·d+
· d = η/k.
2dk
2dk
Finally, we have from upper bound of the error in the individual blocks that
X
f′ (x)
^
∇xi Fj (x) − ∇
F ′ (x) − F
≤
max
xi Fj (x) ∞→∞ ≤ η.
∞→∞
j∈[k]

5

i∈[k]

√
Sample Complexity Upper Bounds with Ω( log k) Separation

√
We now show that a mean separation of order Ω( log k) suffices to learn the model parameters up
to arbitrary accuracy δ > 0, with poly(d, k, log(1/δ)) samples. In all the bounds that follow, the
interesting settings of parameters are when ρ, 1/wmin ≤ poly(k). We note that these upper bounds
are interesting even in the case of uniform mixtures: wi , σi being equal across the components.
For sake of exposition, we will restrict our attention to the case when the standard deviations
σi , and weights wi are known for all i ∈ [k]. We believe that similar techniques can also see used
to handle unknown σi , wi as well (see Remark 4.6). In what follows ρσ corresponds to the aspect
ratio of the covariances i.e., ρσ = maxi∈[k] σi / mini∈[k] σi .
Theorem 5.1 (Same as Theorem 1.3). There exists a universal constant c > 0 such that suppose
we are given samples from a mixture of spherical Gaussians G = {(wi , µi , σi ) : i ∈ [k]} (with known
weights and variances) that are ρ-bounded and the means are well-separated i.e.
p
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ c log(ρσ /wmin )(σi + σj ),
(49)

there is an algorithm that for any δ > 0, uses poly(k, d, ρ, log(1/wmin ), 1/δ) samples and recovers with
up to δ relative error i.e., finds {µ′i : i ∈ [k]} such that
 high probability the means

∆param G, {(wi , µ′i , σi ) : i ∈ [k]} ≤ δ.

Such results are commonly referred to as polynomial identifiability or robust identifiability results.
We can again assume as in Section 4 that without loss of generality that d ≤ k due to the following
dimension-reduction technique using PCA [45]. Theorem 5.1 follows in a straightforward manner
by combining the iterative algorithm, with initializers given by the following theorem.
31

Theorem 5.2 (Initializers Using Polynomial Samples). For any constant c ≥ 10, suppose we are
given samples from a mixture of spherical Gaussians G = {(wi , µi , σi ) : i ∈ [k]} that are ρ-bounded
and the means are well-separated i.e.
p
(50)
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ 4c log(ρσ /wmin )(σi + σj ).

There is an algorithm that uses poly(kc , d, ρ) samples and with high probability learns the parameters
of G up to k−c accuracy, i.e., finds another mixture of spherical Gaussians G̃ that has parameter
distance ∆param (G, G̃) ≤ k−c .

The key difference between Theorem 5.1 and Theorem 5.2 is that in the former, the parameter
estimation accuracy is independent of the separation (the constant c > 0 does not
√ depend on δ). If
ρs = kO(1) and wmin ≥ k−O(1) , then we need means µj , µℓ to be separated by Ω( log k)(σj + σℓ ) to
get reasonable estimates of the parameters with poly(k) samples. While the algorithm is sample
2
efficient, it takes time that is (ρ/wmin )O(ck ) time since it runs over all possible settings of the O(k)
parameters in d ≤ k dimensions. Note that the above theorem holds even with unknown variances
and weights that are unequal.
We can use Theorem 5.2 as a black-box to obtain a mixture G (0) such that ∆param (G, G (0) ) ≤ k−c
whose parameters will serve as initializers for the iterative algorithm in Section 4. Theorem 5.2
follows by exhibiting a lower bound on the statistical distance between any two sufficiently separated
spherical Gaussian mixtures which have non-negligible parameter distance.
Proposition 5.3. Consider any two spherical Gaussian mixtures G, G ∗ in Rd as in Theorem 5.2,
with their corresponding p.d.fs being f and f ∗ respectively. There is a universal constant c′ > 0,
1
, and G is well-separated:
such that for any c ≥ 5, if the parameter distance ∆param (G, G ∗ ) ≥ kc−1
p
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ 4c log(ρs /wmin )(σi + σj ),

and both mixtures have minimum weight wmin ≥ 1/kc , then
Z
|f (x) − f ∗ (x)| dx ≥
kf − f ∗ k1 =
Rd

where ρs = max

n max

σj∗ maxj∈[k] σj
minj∈[k] σj , minj∈[k] σj∗
j∈[k]

c′
dk2c ρ2s

,

(51)

o
.

In the above proposition, ρs ≤ ρ2 is a simple upper bound since we will only search for all
parameters of magnitude at most ρ. But it can be much smaller if we have a better knowledge
of the range of {σj : j ∈ [k]}. We note that in very recent independent work, Diakonikolas et
al. established a similar statement about mixtures of Gaussians where the components have small
overlap (see Appendix B in [20]). We first see how the above proposition implies Theorem 5.2.

5.1

Proposition 5.3 to Theorem 5.2

The following simple lemma gives a sample-efficient algorithm to find a distribution from a net
of distributions T that is close to the given distribution. This tournament-based argument is a
commonly used technique in learning distributions.

32

Lemma 5.4. Suppose T is a set of probability distributions over X , and we are given m samples
from a distribution D, which is δ close to some distribution D ′′ ∈ T in statistical distance, i.e.,
kD − D ′′ kT V ≤ δ. Then there is an algorithm that uses m = O(δ−2 log|T |) samples from D and
with probability at least 1 − 1/|T | finds a distribution D ∗ such that kD − D ∗ kT V ≤ 4δ.
Proof. Let T = {D1 , D2 , . . . , D|T | } be the set of distributions over X . For any i 6= j let Aij ⊂ X
be such that
P [x ∈ Aij ] − P [x ∈ Aij ] = kDi − Dj kT V .
(52)
x←Di

x←Dj

The algorithm is as follows.
• Use m = O(δ−2 log|T |) samples to obtain estimates pij satisfying with probability at least
1 − 1/|T | that
∀i, j, |pij − P [x ∈ Aij ]| ≤ δ/2 .
(53)
x←D

• Output the first distribution Di ∈ T that satisfies
∀j ∈ {1, . . . , |T |}, pij − P [x ∈ Aij ] ≤ 3δ/2 .
x←Di

(54)

First, notice that by (53) and the assumption that kD−D ′′ kT V ≤ δ, D ′′ satisfies the test in (54).
We next observe by the definition of Aij in (52) that for any i 6= j such that both Di and Dj pass
the test, we must have kDi − Dj kT V ≤ 3δ. This implies that the output of the algorithm must be
within statistical distance 4δ of D, as desired.
We now prove Proposition 5.2 assuming the above Proposition 5.3. This follows in a straightforward manner by using the algorithm from Lemma 5.4 where T is chosen to be a net over all
possible configuration of means, variances and weights. We give the proof below for completeness.
Proof of Proposition 5.2. Set δ := c′ /(8ρ8 dkc+1 ), where c′ is given in Proposition 5.3, and ε =
δ/(6kdρ). We first pick T to be distributions given by an ε-net over the parameter space, and the
algorithm will only output one of the mixtures in the net. Each Gaussian mixture has a standard
deviation σj ∈ R and k components each of which has a mean in Rd , and weight in [0, 1]. Further,
all the means, variances are ρ bounded. Hence, considering a net T of size M ≤ (ρ/(wmin ε))(d+2)k of
ρ-bounded well-separated mixtures with weights at least wmin , we can ensure that every Gaussian
mixture is ε close to some ρ bounded mixture in the net in parameter
distance ∆param . This
√
corresponds to a total variation distance distance of at most δ = 6k dρε as well, using Lemma A.3.
′′
Hence, there is a mixture
√ of well-separated ρ-bounded spherical Gaussian mixture G ∈ T such
′′
that kG − GkT V ≤ 6k dρε = δ.
By using Lemma 5.4, we can use m = O(δ−2 log|T |) = O(k2c+3 (d + 2)3 ρ16 log(ρkd/wmin )) and
with high probability, find some well-separated ρ-bounded mixture of spherical Gaussians G ∗ such
that kG − G ∗ kT V ≤ 4δ. Proposition 5.3 implies that the parameters of G ∗ satisfy ∆param (G, G ∗ ) ≤
k−c .

33

5.2

Proof of Proposition 5.3

To show Proposition 5.3, we will consider any two mixtures of well-separated Gaussians, and show
that the statistical distance is at least inverse polynomial in k. This argument becomes particularly
tricky when the different components can have different values of σi e.g., instances where one
component of G with large σi is covered by multiple components from G ∗ with small σj∗ values.
For each component j ∈ [k] in G, we define the region Sj around µj , where we hope to show a
statistical distance.
Definition 5.5 (Region Sj ). In the notation above, for the component of G centered around µj ,
define ebjℓ as the unit vector along µ∗ℓ − µj .
p
Sj = {x ∈ Rd : ∀ℓ ∈ [k] |hx − µj , ebjℓ i| ≤ 2c log(ρσ /wmin )σj }.
(55)
The following lemma shows that most of the probability mass from jth component around µj
is confined to Sj .

Lemma 5.6.
∀j ∈ [k],

Z

Sj

gµj ,σ (x)dx ≥ 1 −



wmin
ρσ

4c

.

(56)

Proof. The region Sj is defined by k equations, one for each of the components ℓ in G ∗ . When x is
drawn from the spherical Gaussian gµj ,σj
P

x←Gµj ,σj

h

|hx − µj , ebjℓ i| > 2cσj

p

 2



i
1 wmin 4c
wmin 2c
,
≤
log(ρσ /wmin ) ≤ Φ0,1 (2c log(ρσ /wmin )) ≤
ρσ
k
ρσ

where the last step follows since c ≥ 5 and wmin ≤ 1/k. Hence, performing a union bound over the
k components in G ∗ completes the proof.
The following lemma shows that components of G ∗ that are far from µj do not contribute much
ℓ1 mass to Sj .
Lemma 5.7. For
in G, and let ℓ ∈ [k] be a component of G ∗ such that
p a component j ∈ [k]
∗
∗
kµj − µℓ k2 ≥ 2c log(ρσ /wmin )(σj + σℓ ). Then the total probability mass from the ℓth component
of G ∗ is negligible, i.e.,
2c2

Z
∗ wmin
∗ ∗ ∗
.
(57)
wℓ gµℓ ,σℓ (x)dx < wℓ
ρσ
Sj
Proof. Let ebjℓ be the unit vector along µ∗ℓ − µj . We will now show that every point x ∈ Sj , is far
from µ∗ℓ along the direction ejl :
p
x ∈ Sj =⇒ |hx − µj , ebjℓ i| ≤ 2c log(ρσ /wmin ) · σj
p
p
kµj − µ∗ℓ k2 ≥ 2c log(ρσ /wmin )(σj + σℓ∗ ) =⇒ hµ∗ℓ − µj , ebjℓ i ≥ 2c log(ρσ /wmin )(σj + σℓ∗ )
p
Hence ∀x ∈ Sj hµ∗ℓ − x, ebjℓ i ≥ 2c log(ρσ /wmin ) · σℓ∗
34

Hence the ℓ1 contribution from gµ∗ℓ ,σ∗ restricted to Sj is bounded as
Z
Z
∗ ∗ ∗
wℓ∗ gµ∗ℓ ,σℓ∗ (x) dx
wℓ gµℓ ,σℓ (x)dx ≤
√
Sj
ejℓ i≥2c log(ρσ /wmin )σℓ∗
x:hµ∗ℓ −x,b

2c2
∗ wmin
∗
≤ wℓ Φ0,1 (2c log(ρσ /wmin )) ≤ wℓ
,
ρσ
as required.
The following lemma shows that there is at most one component of G ∗ that is close to the
component (wj , µj , σj2 I) in G.
Lemma 5.8 (Mapping between centers). Suppose we are given two spherical mixtures of Gaussians
G, G ∗ as in Proposition 5.3. Then for every j ∈ [k] there is at most one ℓ ∈ [k] such that
p
(58)
kµj − µ∗ℓ k2 ≤ 4c log(ρσ /wmin )σℓ∗ .

Proof. This follows from triangle inequality. Suppose for contradiction that there are two centers
corresponding to indices
p ℓr with r = 1, 2, that satisfy (58). By using triangle inequality this shows
that kµ∗ℓ1 − µ∗ℓ2 k ≤ 4c log(ρσ /wmin )(σℓ∗1 + σℓ∗2 ) which contradicts the assumption.

We now proceed to the proof of the main proposition (Proposition 5.3) of this section. We will
try to match up components in G, G ∗ that are very close to each other in parameter distance and
remove them from their respective mixtures. Then we will consider among unmatched components
the one with the smallest variance. Suppose (wj , µj , σj ) were this component, we will show a
significant statistical distance in the region Sj around µj .
The following lemma considers two components, Gj = (wj , µj , σj ) from G , and G∗j = (wj∗ , µ∗j , σj∗ )
from G ∗ that have a non-negligible difference in parameters (we use the same index j for convenience,
since this is without loss of generality). This lemma shows that if σj ≤ σj∗ , then there is some region
S ⊂ Sj where the component Gj has significantly larger probability mass than G∗j . We note that
it is crucial for our purposes that we obtain non-negligible statistical distance in a region around
Sj . Suppose fj and fj∗ are the p.d.f.s of the two components (with weights), it is easier to lower
bound kfj − fj∗ k1 (e.g. Lemma 38 in [35]). However, this does not translate to a corresponding
lower bound restricted to region S i.e., kfj − fj∗ k1,S since fj and fj∗ do not represent distributions
(e.g. kfj k1 = wj < 1).
Lemma 5.9. For some universal constant c1 > 0, suppose we are given two spherical Gaussian
components with parameters (wj , µj , σj ) and (wj∗ , µ∗j , σj∗ ) that are ρσ -bounded satisfying
σj ≤ σj∗ and

kµj − µ∗j k2 |σj − σj∗ |
+
+ |wj − wj∗ | ≥ γ.
σj
σj

Then, there exists a set S ⊂ Sj such that
Z
c1 γ 2
.
wj gµj ,σj (x) − wj∗ gµ∗j ,σj∗ (x) dx >
dρ2s
S

(59)

(60)

Before we proceed, we present two lemmas which lower bound the statistical distance when the
means of the components differ, or if the means are identical but the variances differ.
35

Figure 1: Regions Sj used in Lemma 5.10 and Lemma 5.11: Lemma 5.10 shows a statistical difference
in one of the strips Lj or Rj . Lemma 5.11 shows a statistical difference in one of the annular strips R1 or R2
Lemma 5.10 (Statistical Distance from Separation of Means). In the notation of Lemma 5.9,
suppose kµj − µ∗j k2 ≥ γσj . Then, there exists a set S ⊂ Sj such that
Z
γσj2
γ
∗ ∗ ∗
(61)
wj gµj ,σj (x) − wj gµj ,σj (x) dx >
2 > 8ρ2 .
∗
8(σj )
S
s
Proof. Without loss of generality we can assume that µj = 0 (by shifting the origin to µj ). Let us
consider two regions
Rj = Sj ∩ {x : hx, ebjj i ∈ [ √12π σj , √22π σj ]} and Lj = Sj ∩ {x : hx, ebjj i ∈ [− √22π σj −

Due to the symmetry of Sj (Definition 5.5), it is easy to see that
Z
Z
1
f (x)dx =
f (x)dx ≥ .
8
Lj
Rj
On the other hand, we will now show that
Z
f ∗ (x)dx >

1+

2γσj2

Rj

σj∗ 2

!Z

f ∗ (x)dx

√1 σj ]}.
2π

(62)

(63)

Lj

For any point x, let xjl = hx, ebjj i. Further, µ∗j = kµ∗j kb
ejj . We first note that x ∈ Rj ⇐⇒
(−x) ∈ Lj due to the symmetric definition of Sj . Hence,
!
Z
Z
Z
∗ (−x)
f
−π
f ∗ (x)dx =
f ∗ (x)dx =
exp
(k−x − µ∗j k22 − kx − µ∗j k22 ) f ∗ (x)dx
∗ (x)
∗2
f
σ
Lj
Rj
Rj
j
Z
Z




∗
∗2
∗
exp −2πkµj k2 hx, ebjj i/(σj∗ )2 f ∗ (x)dx
exp −2πhx, µj i/σj f (x)dx =
=
Rj
Rj
Z
√
σj
f ∗ (x)dx, since hx, ebjj i ≥ √
≤ exp(− 2πσj2 γ/σj∗ 2 )
for x ∈ Rj .
2π
Rj
36

Hence from (62) and (63) either
R
R
√
√
∗
2πγσj2
2πγσj2
Rj f (x)dx
Lj f (x)dx
R
R
or
,
>
1
+
>
1
+
∗
σj∗ 2
σj∗ 2
Lj f (x)dx
Rj f (x)dx
R
R
which completes the proof, since Lj f (x)dx = Rj f (x)dx ≥ 1/8.

Lemma 5.11 (Statistical Distance from Different Variances). In the notation of Lemma 5.9, suppose µj = µ∗j , but σj = (1 − η)σj∗ . Then, there exists a set S ⊂ Sj such that
Z
wj gµj ,σj (x) − wj∗ gµ∗j ,σj∗ (x) dx > min{c2 η, 1},
(64)
S

where c2 > 0 is a universal constant.
Proof. Without loss of generality, let us assume that µj = µ∗j = 0 (by shifting the origin), and
√
√
√
√
σj = 1 (by scaling). Let R1 = [ √12π ( d − 2), √12π ( d − 1)] and R2 = [ √12π ( d + 1), √12π ( d + 2)]
and consider two annular strips T1 = {x : kxk ∈ R1 } and T2 = {x : kxk ∈ R2 }.
First, we note that using standard facts about the χ2 (d) distribution,for some appropriately
chosen universal constant c2 > 0 we have
Z
Z
f (x)dx ≥ c2 .
f (x)dx ≥ c2 , and
T2

T1

We will show that there is a significant statistical difference between f, f ∗ in either T1 or T2 .
Let us assume for contradiction that
R
R
f (x)dx
T1 f (x)dx
R T2
≤
(1
+
η)
and
(1
−
η)
≤
≤ (1 + η).
(65)
(1 − η) ≤ R
∗
∗
T1 f (x)dx
T2 f (x)dx
Let us consider the range of values that f takes in T1 and T2 . We


 √
 √
wj
wj
∀x ∈ T1 d exp −( d − 1)2 /2 ≤ f (x) ≤ d exp −( d − 2)2 /2
σj
σj
 √
 √


wj
w
j
∀x ∈ T2 d exp −( d + 2)2 /2 ≤ f (x) ≤ d exp −( d + 1)2 /2
σj
σj



√
√
√
f (x1 )
≥ exp − 12 ( d − 1)2 ) − ( d + 1)2
≥ exp(2 d) ≥ e2 .
∀x1 ∈ T1 , x2 ∈ T2 ,
f (x2 )

Let A(r) be the d − 1 dimensional volume of {x ∈ Sj : kxk = r}. From (65), we have that
R
R
wj r∈R1 exp(−r 2 /2)A(r) dr
T1 f (x)dx
R
R
= ∗
d
2 2
∗ (x)dx
w
(1
−
η)
f
j
r∈R1 exp(−(1 − η) r /2)A(r) dr
T1
R
wj r∈R1 exp(−r 2 /2)A(r) dr
R
= ∗
wj (1 − η)d r∈R1 exp ((2η − η 2 )r 2 /2) · exp(−r 2 /2)A(r) dr


√
wj
2
≥ ∗
·
exp
−η(
d
−
1)
/2
.
wj (1 − η)d
 √

wj
2
≤
(1
+
η)
exp
η(
Hence,
d
−
1)
/2
.
wj∗ (1 − η)d
37

(66)
(67)
(68)

Using a similar argument for T2 we get
 √

wj
2
d
+
1)
≥
(1
−
η)
exp
η(
/2
.
wj∗ (1 − η)d
Combining the two we get

 √

 √
2
2
(1 + η) exp η( d − 1) /2 ≥ (1 − η) exp η( d + 1) /2
 √ 
exp 2η d ≤ (1 + η)2 ≤ e2η ,

which is a contradiction.

Proof of Lemma 5.9. Let fj (x) = wj gµj ,σj (x), fj∗ (x) = wj∗ gµ∗j ,σj∗ (x) be the p.d.f. of the two components. From (59) we know that there is some non-negligible separation in the parameters. Hence
either the weights, or means, or variances are separated by at least Ω(γ). We will now consider
three cases depending on whether there is non-negligible separation in the means, variances or
weights (in that order).
Set γ1 := γ 2 /(256d), γ2 := c2 γ1 /2 where c2 is the constant in Lemma
5.11. Suppose there is
√
some non-negligible separation in the means i.e., kµj − µ∗j k ≥ γ2 σj / 2π. Lemma 5.10 shows that
there is a set S ⊂ Sj that has
γ2 σj2
c1 γ 2
,
kf − f ∗ k1,S ≥
≥
dρ2s
8σj∗ 2
√
where c1 = c2 /(512 2π).
√
Otherwise we have that kµj − µ∗j k < γ2 σj / 2π. Let gj∗ be the p.d.f. of the Gaussian component
(wj∗ , µj , σj∗ ). From Lemma A.3 we know that kfj∗ − gj∗ k1 ≤ γ2 . Now, fj and gj∗ are p.d.f. of
components that have the same mean.
If σj∗ − σj > γ1 σj . From Lemma 5.11 we see that
kf − g ∗ k1,S ≥ c′ γ1 =⇒ kf − f ∗ k1,S ≥ c′ γ1 − γ2 ≥ c′ γ1 /2.
Finally if σj∗ − σj ≤ γ1 σj , then one can bound the statistical distance over S = Sj between two
components with equal means and variances, but different weights:
Z
wj gµj ,σj (x) − wj∗ gµj ,σℓ (x) dx − γ2
kfj − fj∗ k1,S ≥ kfj − gj∗ k1,S − kgj∗ − fj∗ k =
S
Z
Z
∗
wj∗ gµj ,σj (x) − wj∗ gµj ,σj∗ (x) dx − γ2
wj gµj ,σj (x) − wj gµj ,σj (x) dx −
≥
S
S
Z
p
∗
≥ |wj − wj | gµj ,σj (x) − 2 dγ1 − γ2
( from Lemma A.3)
S

 γ γ
2γ
1 − (wmin /ρσ )4c − − ≥ γ/4,
≥
3
8
8

where the last line follows from Lemma 5.6.

We now complete the proof of the main proposition of this section, that lower bounds the statistical difference between two mixtures of well-separated Gaussians which differ in their parameters.
38

Proof of Proposition 5.3. We follow the approach outlined earlier in the section. We first match
up components from (wj , µj , σj ) ∈ G and (wj∗ , µ∗j , σj∗ ) ∈ G ∗ that are at most γ = k−c close in
parameter distance. From triangle inequality and well-separatedness of G, G ∗ , it is easy to see that
one component from G (or G ∗ ) cannot be matched to more than one component from G ∗ (or G
respectively). Let {1, 2, . . . , k ′ } be the indices of the components in G and G ∗ that are unmatched,
for ℓ ∈ {k′ + 1, . . . , k}, let the component (wℓ , µℓ , σℓ ) of G be matched to the component (wℓ∗ , µ∗ℓ , σℓ∗ )
of G ∗ . Since the parameter distance ∆param (G, G ∗ ) ≥ k−c+1 , we have that k′ ≥ 1.
Consider among the unmatched components in both G and G ∗ , the one with the smallest
variance: let this component be (wj , µj , σj ) from G without loss of generality. From Lemma 5.8,
we know that at most one component of G ∗ satisfies (58). Again, without loss of generality, let
(wj∗ , µ∗j , σj∗ ) be this component of G ∗ (if it exists). Hence
p
p
∀ℓ ∈ [k′ ], ℓ 6= j, kµj − µ∗ℓ k2 > 4c log(ρσ /wmin )σℓ∗ ≥ 2c log(ρσ /wmin )(σj + σℓ∗ ),

where the last inequality follows because σj2 is the smallest variance.4 Further, all the matched
components ℓ ∈ {k′ + 1, . . . k} in G ∗ are all far from µj since
∀ℓ ∈ {k′ + 1, . . . , k}, kµj − µ∗ℓ k2 ≥ kµj − µℓ k2 − kµℓ − µ∗ℓ k
p
≥ 4c log(ρσ /wmin )(σj + σℓ ) − γ
p
p
≥ 4c log(ρσ /wmin )(σj + σℓ∗ ) − γ(1 + 4c log(ρσ /wmin ))
p
≥ 2c log(ρσ /wmin )(σj + σℓ ),

since γ < 2/(5ρσ ). From Lemma 5.7, there is negligible contribution from the rest of the components
(using c ≥ 3):


Z
X
wmin 4c
∗
∗
∗
.
(69)
gµℓ ,σℓ (x)dx <
wℓ
ρσ
Sj
ℓ∈[k]:ℓ6=j

From Lemma 5.9 there is a subset S ⊂ Sj where there is significant statistical distance
Z

S

wj gµj ,σj (x) − wj∗ gµ∗j ,σj∗ (x) dx >

c1 γ 2
.
dρ2s

Combining the last two equations, we have
∗

kf − f k1 ≥
≥
≥

Z

S

Z

S

X
ℓ

wℓ gµℓ ,σℓ (x) −

wj gµj ,σj (x) −

c1 γ 2
−
dρ4s



1
ρσ k

4c

X

wℓ∗ gµ∗ℓ ,σℓ∗ (x) dx

ℓ

wj∗ gµ∗j ,σj∗ (x)
≥

c′ γ 2
dρ2s

dx −

XZ
ℓ6=j

S

wℓ∗ gµ∗ℓ ,σℓ∗ (x)dx

for some universal constant c′ > 0, since γ ≥ k−c .
4

This is in fact the only point where we use the fact that we picked the small variance Gaussian.

39

6

Efficient Algorithms in Low Dimensions

In this section, we give a computationally efficient algorithm that works in d = O(1) dimensions,
and learns the mixture of k spherical Gaussians even when the separation
pbetween centers is O(σ).
In comparison, previous algorithms need separation of the order of Ω(σ log(k/δ)). We prove the
following theorem.
Theorem 6.1. There exists universal constants c > 0 such that the following holds. Suppose we
are given samples from a mixture of spherical Gaussians G = {(wj , µj , σj ) : j ∈ [k]}, where the
weights and covariances are known, such that kµj k ≤ ρ ∀j ∈ [k] and
√

p
∀i, j ∈ [k], i 6= j : kµi − µj k2 ≥ c
d + log(ρσ ρw ) · (σi + σj ).
(70)

O(d)
−1
For any δ > 0, there is an algorithm using time (and samples) poly wmin
, δ−1 , ρ, ρσ
that
with high probability recovers the means up to δ accuracy i.e. finds for each j ∈ [k], µ̃j such that
kµ̃j − µj k2 ≤ δσj .

In the above theorem, when both √
ρw , ρσ = O(1) as in the case of uniform mixtures, this
corresponds to a separation of order Ω( d).
The above theorem follows by applying the guarantees of the iterative algorithm (Theorem 4.1)
along with a computationally efficient procedure that finds appropriate initializers. The following
theorem shows how to find reasonable initializers for µj , σj , wj for each of the k components.
Theorem 6.2. Let c0 ≥ 2 be any constant, and suppose ε0 = exp(−c0 d). There is an algorithm
 O(d)
running in time ρd
poly(1/wmin ) that given samples from a ρ-bounded mixture of k spherical
ε30
Gaussians G = {(wj , µj , σj ) : j ∈ [k]} in d dimensions satisfying

√
p
d + log(ρw ρσ ) (σi + σj ),
(71)
∀i 6= j ∈ [k], kµi − µj k2 ≥ 4c0
can find with high probability {(µ̃j , σ̃j , w̃j ) : j ∈ [k]} s.t.
√
∀j ∈ [k], kµ̃j − µj k2 ≤ ε0 σj d, |σ̃j − σj | ≤ ε0 σj , and |w̃j − wj | ≤ ε0 wj .

Note that the above theorem also finds initializers for the weights and variances when they are
unknown. Hence, Theorem 6.1 will also apply to the setting with unknown weights and variances
if we get similar guarantees for Theorem 4.1 (see Remark 4.6).
To prove Theorem 6.2, we will first find estimates µ̃j for the means µj (Proposition 6.3), and
then obtain good estimates σ̃j , w̃j for each j ∈ [k]. Proposition 6.3 already suffices in the case of
known weights and variances.
O(d)

Proposition 6.3. In the notation of Theorem 6.2, there is an algorithm running in ε3 wρ
0 min
√
time that finds w.h.p. µ̃1 , . . . , µ̃k s.t. ∀j ∈ [k], kµ̃j − µj k2 ≤ ε0 σj d.

40

We first start with expressions for the first derivative (gradient) and second derivative (Hessian)
at a point x in terms of the model parameters.
′

f (x) = ∇f (x) = −2π
f ′′ (x) = ∇2 f (x) = 4π 2

k
X
wj
j=1

σj

k
X
wj
j=1

σj2

· gµj ,σj (x) ·
· gµj ,σj (x) ·

(x − µj )
σj
!
1
1
Id×d .
(x − µj )⊗2 −
2π
σj2

(72)

(73)

Note that using polyd (k, ρ) samples, we have access to the p.d.f. f (x) and its derivatives f ′ = ∇f ,
f ′′ = ∇2 f up to 1/polyd (k) accuracy at any point x.
The algorithm will consider a δ-net of points, and find “approximate local-maxima” of the p.d.f.,
which are defined as follows.
Definition 6.4 (Approximate Local Maximum). Consider a mixture of k spherical Gaussians in Rd
with parameters {(wj , µj , σj ) : j ∈ [k]} and p.d.f. f . Then x ∈ Rd is an approximate local-maxima
iff
√
πε0 dσmin
wmin
′
(ii) kf (x)k2 ≤
f (x),
(74)
(i) f (x) ≥ d ,
2
2σmax
4σmax
π
(iii) f ′′ (x) = ∇2 f (x)  − 2 f (x)Id×d
.
2σmax
To show the above proposition, we will show that all approximate local-maxima are close to
one of the means µj (Lemma 6.5 and Lemma 6.6), and there is at least one such approximate local
maxima near each mean µj (Lemma 6.7). Further, since the parameters are separated, this will
allow us to pick all approximate local-maxima in a net, cluster them geometrically and pick one
such point from each cluster to get good initializers for each mean µj . These statements will also
allow for some slack to tolerate estimation errors.
We start with some inequalities that use the separation between means to show the p.d.f. and
the first few “moments” near one of the means µj is dominated by the the jthpcomponent. By
applying Lemma 4.13 we get that at any point x ∈ Rd such that kx − µj k2 ≤ σj d/π

σ
X
min 2
wi gµi ,σi (x) < wj gµj ,σj (x) exp(−c0 d)
.
(75)
σmax
i6=j∈[k]

X

i6=j∈[k]

X

i6=j∈[k]

wi ·
wi ·

σ

kx − µi k2
min 4
· gµi ,σi (x) < wj gµj ,σj (x) exp(−c0 d) ·
.
σi
σmax


σ
kx − µi k22
min 4
.
·
g
(x)
<
w
g
(x)
exp(−c
d)
·
µi ,σi
j µj ,σj
0
σmax
σi2

(76)
(77)

The following lemma shows that any point that is far from all of the means is not an approximate
local maximum (does not satisfy condition (iii) of Def. 6.4).
Lemma 6.5. In the notation of Proposition 6.3, for any x ∈ Rd that satisfies
r
d
, ∃u ∈ Rd s.t. uT Hx u > 0,
∀j ∈ [k] kx − µj k2 > σj
π
41

where where Hx = f ′′ (x) represents the Hessian evaluated at x. Hence, such a point x is not an
approximate local maxima.
Proof. Let Hx,j = (x − µj )⊗2 − σj2 I/(2π). Then
dσj2
dσj2
>
∀j ∈ [k], tr(Hx,j ) = kx − µj k −
2π
2π
2

kx − µj k2
since
>
σj

r

d
.
π

Let v a random unit standard Gaussian vector drawn from N (0, 1)d .

!
k
E[v T Hx v]
4π 2 X wj
tr(Hx )
πkx − µj k2
=
=
· exp −
· tr(Hj )
d
d
E[kvk22 ]
σj2
σ d+4
j=1 j
!
k
πkx − µj k2
d
1
2π X wj
exp −
· 2 > 2π · f (x) ·
>
2
2.
d
d
σ
σ
max
σ
j σj
j
j
j=1 j

2 .
Hence, there is a direction v T Hx v ≥ 2πf (x)kvk22 /σmax

The following lemma shows that we cannot have approximate
plocal maxima (or more generally,
√
critical points) whose distance from µj is between [ε0 dσmin , d/π · σj ].√ Hence, together with
Lemma 6.5, this shows that every approximate local maximum is within ε0 dσmin from one of the
true means.
Lemma 6.6. In the notation of Proposition 6.3,
p 0 d) with c0 ≥ 2, suppose
√ for any ε0 = 2 exp(−c
x ∈ Rd satisfies for some component j ∈ [k], ε0 dσmin < kx − µj k2 ≤ d/πσj , then
√
√
ε0 dσmin
ε0 dσmin
′
kf (x)k2 = k∇f (x)k2 >
· f (x).
· f (x) ≥
2
4σmax
4σj2
Proof. From (72), the first derivative satisfies
kf ′ (x)k2 ≥

X 2πwi kx − µi k2
2πwj kx − µj k2
·
· gµi ,σi (x).
· gµj ,σj (x) −
2
σi
σi
σj

(78)

i6=j∈[k]

Applying (76) with the given separation,

σ
X wi kx − µi k2
ε0 wj σmin
1
min 2
<
·
·gµi ,σi (x) <
·wj gµj ,σj (x) exp(−c0 d)·
gµj ,σj (x). (79)
2
σi
σi
σmin
σmax
2σmax

i6=j∈[k]

√
Further, kx − µj k2 ≥ ε0 dσj . Using (78) and (79),

k

kf ′ (x)k2 ≥ 2πwj gµj ,σj (x) ·

kx − µj k22 X
kx − µi k22
−
2πwi gµi ,σi (x) ·
2
σj
σi2
i6=j

2πkx − µj k2
πε0 wj σmin
≥
· wj gµj ,σj (x) −
gµj ,σj (x) (from (79)),
2
2
σmax
σj
√
πε0 dσmin
· f (x),
≥
2
4σmax
√
where the last inequality follows from (75) and using kx − µj k2 > ε0 dσmin .
42

We now proceed to the proof of Lemma 6.7, which shows that any point that is sufficiently close
to one of the component√means is an approximate local maxima. This shows that in any δ-net (for
2 /σ
sufficiently small δ < ε0 dσmin
max ) , there will be an approximate local maxima.
Lemma 6.7. In the notation of Proposition 6.3, for any ε′ ≤ exp(−c0 d) with c0 ≥ 2 and ∀j ∈ [k],

any point x ∈ Rd with kx − µj k2 ≤

√
ε′ d
32

·

σj2
σmax

is an approximate local maxima. In particular,
√
4πε′ d
(ii) kf (x)k2 ≤
f (x),
σmax

3wmin
(i) f (x) ≥
,
d
4σmax

′

(iii) f ′′ (x) = ∇2 f (x)  −

3π
π
f (x)I  − 2 f (x)I
4σmax
σj2

(80)

.

Proof of Lemma 6.7. The lemma follows in a straightforward way from (75), (76), (77), since
f, f ′ , f ′′ at x are dominated by the jth component. Firstly by considering just the contribution to
the p.d.f. from the jth component, the lower bound on f (x) follows. Now we bound kf ′ (x)k2 .
k
X

kx − µi k2
kx − µj k2
≤ 4πwj gµj ,σj (x) ·
(from (76)).
2
σi
σj2
i=1
√
√
4πε′ d
4πε′ d
≤
· wj gµj ,σj (x) ≤
· f (x).
σmax
σmax

kf ′ (x)k2 ≤ 2π

wi gµi ,σi (x) ·

We argue about f ′′ (x) similarly. Suppose we use M to denote the following matrix, and kM k to
represent its maximum singular value,


X
1
1
1
⊗2
2
(x − µi ) −
Id×d .
wi gµi ,σi (x) · 2
M = 4π
2π
σi σi2
i6=j


X
1
1 kx − µi k22
2
+
kM k ≤ 4π
wi gµi ,σi (x) · 2
2π
σi
σi2
≤

i6=j
2
4π

2
σmax

wj gµj ,σj (x) exp(−c0 d) <

π
wj gµj ,σj (x),
2σj2

where the last line follows from (75), (77) and since exp(−c0 d) < 1/(8π). Further 2πkx−µj k22 /σj2 <
1/4. Substituting in (73), we get
!
2
2πkx
−
µ
k
2π
j 2
f ′′ (x)  2 wj gµj ,σj (x) −I +
I + kM kI
σj
σj2
 −π
−3π
2π
f (x),
 − 2 wj gµj ,σj (x) I − 14 I − 14 I  2 wj gµj ,σj (x) 
σj
σj
4σj2
where the last inequality follows from (75).
We now proceed to the algorithm and proof of Proposition 6.3.
Proof of Proposition 6.3. The algorithm first considers a δ-net Xδ in Rd over a ball of radius 2ρ, and
√ 3
−(d+4)
2
estimates f (y) up to additive accuracy γ where γ = wmin σmax ε30 δ2 /4 and δ = ε0 dσmin
/(64σmax
)
43

will suffice5 . Similarly, we can also estimate f ′ (y) in ℓ2 norm, and f ′′ (y) in operator norm within
additive accuracy γ. The size of the net is |Xδ | = (ρ/δ)O(d) , and the sample complexity is O(1/γ 2 ) ·
(ρ/δ)O(d) samples.
From Lemma 6.6 and Lemma 6.5 we have that if
√
kf ′ (x)k2
πε0 dσmin
π
−d
′′
≤
,
(81)
f (x) ≥ wmin σmax /2, f (x)  − 2 f (x)I, and
2
2σmax
f (x)
4σmax
√
then there exists j ∈ [k] s.t. kx − µj k2 ≤ ε0 dσmin . √
On the other hand, applying Lemma 6.7 with
3 /σ 2 ) close to µ satisfy
ε′ = ε0 σmin /(32σmax ), any point that is within O(ε0 dσmin
j
max
√
πε0 dσmin
3π
kf ′ (x)k2
3wmin
′′
2
≤
.
(82)
f (x) ≥ d , f (x) = ∇ f (x)  − 2 f (x)I, and
2
4σmax
4σmax
f (x)
8σmax
Our accuracy γ of estimating f, f ′ , f ′′ is chosen so that we can distinguish between the bounds
in (81) and (82). For convenience, since we have sufficiently accurate estimates, we will abuse
notation and also use f (x), f ′ (x), f ′′ (x) to represent the estimate of the f, f ′ , f ′′ at x.
First using our estimates, we consider all points
√
n
πε0 dσmin o
kf ′ (y)k2
π
−d
′′
≤
.
T = y ∈ Xδ | f (y) ≥ wmin σmax /2, f (y)  − 2 f (y)I,
2
2σmax
f (y)
4σmax
−(d+2)

We can find T from our estimates since γ < wmin σmax
√ /8. From (81), we have that for every
y ∈ T , there is some j ∈ [k], such that ky − µj k2 ≤ ε0 dσmin .
Further the means are well separated i.e., kµi − µj k2 > 4(σi + σj ). Hence, suppose we define
√
Tj∗ = {y ∈ T | kµj − yk2 ≤ ε0 dσj },
the sets {Tj∗ : j ∈ [k]} are disjoint and form a partition of T that is consistent with the k components
of the Gaussian mixture G. From the separation conditions we see that the distances between any
two points in the same cluster Tj∗ are smaller than the distance between any two points in different
clusters Tj∗1 and Tj∗2 (j1 6= j2 ). Hence, we can use single-linkage clustering algorithm (see Awasthi
et al. [6] for a proof that single-linkage algorithm suffices)to find the clustering {Tj∗ : j ∈ [k]} in
time poly(|T |) ≤ poly(|Xδ |).
Finally, for each j ∈ [k], let µ̃j be any point in Tj∗ . Since, the coarseness of the net δ satisfies
√ 3
2 ), we have from (82) that there is at least one point y ∈ T close to each µ .
/(64σmax
δ < ε0 dσmin
j
√
Hence, kµ̃j − µj k2 ≤ ε0 dσmin , as required.
Note. In fact, the above proposition can also be used to show that f has√exactly k local maxima
r1 , r2 , . . . , rk , such that there is a unique rj satisfying krj − µj k2 ≤ εσj d. This is by using a
quantitative version of the inverse function mapping theorem with the function h(x) = ∇f (x) = 0
(one can use the Newton method as in Section 4).
Again, in what follows, when it is clear that we have sufficiently accurate estimates, we will
abuse notation and also use f (x), f ′ (x), f ′′ (x) to represent the estimate of the f, f ′ , f ′′ at x.
5
For our purposes, it suffices to have estimates of f, f ′ , f ′′ at each of the points of the net Xδ . Hence we can just
take histogram counts in a small ball around each of the net points, and estimate f up to accuracy γ with |Xδ |/δ 2
samples. Similarly the derivatives can also estimated using poly(|Xδ |, 1δ , γ1 ) samples.

44

Lemma 6.8. Assume the conditions in Theorem 6.2, and let ε0 = exp(−c0 d). Suppose we have
O(d)

µ̃j ∈ Rd such that kµ̃j − µj k2 ≤ exp(−2c0 d)σj . Then there is an algorithm running in ε3 wρ
0 min
time, that w.h.p. finds σ̃j ∈ R+ such that |σ̃j − σj | ≤ ε0 σj .

Proof. Let κ be √
a fixed number chosen so that κ ≤ σj , and pick any point y ∈ Rd such that
1
ky − µ̃j k2 = √π κ d.6 Based on the estimates of the p.d.f. at µ̃j , y, we will set
√
κ d
σ̃j = r 
.
f (µ̃j )
log f (y)

√
Both µ̃j , y are in a ball of radius at most σj d/π around µj . Further from Lemma 4.13, and
since we have good estimate for the p.d.f. f (µ̃
 j ) and
 f (y) at these points, we have f (µ̃j ) =
2

wj σj−d (1 ± exp(−2c0 d)), and f (y) = wj σj−d exp − κσ2d
j

log



f (µ̃j )
f (y)



−

(1 ± exp(−2c0 d)). Dividing,

κ2 d
≤ 2 exp(−2c0 d).
σj2

Substituting for σ̃j we have

Hence,

σj2
σ̃j2



f (µ̃ )
σ̃j2 log f (y)j


where |η| ≤ 2 exp(−2c0 d)
σj2 =
f (µ̃ )
log f (y)j + η


η

 .
= σ̃j2 1 +
f (µ̃ )
log f (y)j
− 1 ≤ ε0 , where the last inequality follows from our choice of y, since log(f (µ̃j )/f (y)) ≤

d < exp(c0 d).
Lemma 6.9. Assume the conditions in Theorem 6.2, and let ε0 = exp(−c0 d). Suppose we have
µ̃j ∈ Rd , σ̃j ∈ Rd such that kµ̃j − µj k2 + |σ̃j − σj | ≤ exp(−2c0 d)σ. Then there is an algorithm

O(d)
running in ε3 wρ
time, that w.h.p. finds w̃j ∈ [0, 1] such that |w̃j − wj | ≤ ε0 wj .
0

min

Proof. Let η = wj exp(−2c0 d). Let cd be the constant that is only dependent on d given by
Z

2
cd =
√ exp −πkxk2 dx.
1
x∈Rd :kxk2 ≤ √
2π

d

In fact, cd = γ( d2 − 1, 21 ) is the incomplete Gamma function evaluated at ( d2 − 1, 12 ) which has an
asymptotic approximation given in [21, 8.11(ii)]. Also cd ≥ 2−d/2 /d. To get an estimate of w̃j , we
will consider the set
√
Tj = {x ∈ Rd | kx − µ̃j k2 ≤ √12π dσ̃j }.
6

We can guess such a κ by either doing a binary search, or set it to be one-eighth of the diameter of cluster Tj∗
defined in proof of Proposition 6.3

45

We will now generate N = O(ρ log(dk)/η 2 ) samples x(1) , . . . , x(N ) from the mixture of k Gaussians
and estimate the fraction of samples that are in Tj :7
w̃j =

N
1 X (ℓ)
I[x ∈ Tj ].
cd N

(83)

ℓ=1

From Lemma 4.15, we have small contribution from the other components
|w̃j − E[w̃j ]| = w̃j −

1
cd

Z

y∈Tj

f (y) dy ≤ η ≤ exp(−2cd)wj .

From Lemma 4.13, we have
∀y ∈ Tj , f (y) − wj gµj ,σj (y) ≤ exp(−2c0 d)wj gµj ,σj (y)
Z
Z
wj
wj
wj
Hence, E[w̃j ] −
exp(−2c0 d).
gµj ,σj (y) dy ≤ exp(−2c0 d) ·
gµj ,σj (y) dy <
cd y∈Tj
cd y∈Tj
cd

√
1
dσj }, and let Sd be the surface area of the unit ball in d-dimensions
Let B = {y | ky − µj k2 ≤ 2π
(volume of the d-sphere). Since kµ̃j − µj k + |σj − σ̃j | ≤ exp(−2c0 d)σj , the probability mass on
Tj \ B is small
Z

y∈Tj

gµj ,σj (y) dy −

Z

y∈B

gµj ,σj (y) dy ≤

Z

|

√
2πky−µj k2
−1|≤exp(−2c0 d)
σj

≤ Sd



d
2π

Further, the probability mass inside B is given by

d/2

wj
cd

R

gµj ,σj (y) dy.

× 2 exp(−2c0 d) ≤ exp(−2(c0 − 1)d).

y∈B gµj ,σj (y) dy

= wj . Hence,

Z
Z
Z
wj
wj
wj
gµ ,σ (y) dy +
gµ ,σ (y) dy
gµ ,σ (y) dy −
|E[w̃j ] − wj | ≤ E[w̃j ] −
cd y∈Tj j j
cd y∈B j j
cd y∈Tj j j
wj
wj
exp(−2(c0 − 1)d) +
exp(−2c0 d) ≤ wj exp(−c0 d).
≤
cd
cd

Proof of Theorem 6.2. The proof follows by using Proposition 6.3, followed by Lemma 6.8 and
Lemma 6.9 in that order. Set ε0 = exp(−c0 d). First we use Proposition 6.3 to find w.h.p. initializers
for the means (µ̃j : j ∈ [k]) such that kµ̃j − µj k2 ≤ exp(−4c0 d)σmin . Then using Lemma 6.8, we
find w.h.p. initializers (σ̃j : j ∈ [k]) such that |σ̃j − σj | ≤ exp(−2c0 d)σj . Finally, these initializers
µ̃j , σ̃j can be used in Lemma 6.9 to find w.h.p. w̃j ∀j ∈ [k] such that |w̃j − wj | ≤ exp(−c0 d)wj . By
choosing the failure probability of at most 1/3k2 in each step, we see that the algorithm succeeds
O(d)

.
w.h.p. and runs in time ε3 wρ
0

min

Appendix
7

We could also integrate the estimated p.d.f. over the set Tj to get this estimate.

46

A

Standard Properties of Gaussians

Lemma A.1. Suppose x ∈ R be generated according to N (0, σ 2 ), let Φ̃0,σ (t) represents the probability that x > t, and let Φ̃−1
0,σ (y) represent the quantile t at which Φ̃0,σ (t) ≤ y. Then
t
σ

t2

2

( σt 2 + 1)

e− 2σ2 ≤ Φ̃0,σ (t) ≤

σ − t22
e 2σ .
t

(84)

Further, there exists a universal constant c ∈ (1, 4) such that
p
1p
t
log(1/y) ≤ ≤ c log(1/y).
c
σ

(85)

Lemma A.2. For any σ > 0, q ∈ Z+ and any τ ≥ 2q,
Z
|x|q exp(−πx2 /σ 2 )dx ≤ σ q exp(−2τ 2 ).

(86)

|x|≥τ σ

Lemma A.3. Let p, q correspond to the (weighted) probability density functions of the spherical
Gaussian components in d dimensions with parameters (w1 , µ1 , σ12 ) and (w2 , µ2 , σ22 ) respectively.
Then
!
p
√
r
σ 
2π · kµ1 − µ2 k2
d|σ12 − σ22 |
2
kp − qk1 ≤ |w1 − w2 | + min{w1 , w2 }
. (87)
+
+ 2d ln
σ2
σ2
σ1
Proof. Without loss of generality let w2 ≤ w1 . The KL divergence between any two multivariate
Gaussian distributions with means µ1 , µ2 and covariances Σ1 , Σ2 respectively is given by [22]

 det(Σ ) 
1
2
−1
T −1
.
tr(Σ2 Σ1 ) + (µ1 − µ2 ) Σ2 (µ1 − µ2 ) − d + ln
dKL (N (µ1 , Σ1 )kN (µ2 , Σ2 )) =
2
det(Σ1 )
Applying this to p′ := p/w1 , q ′ := q/w2 we get

 σ 
1 2πkµ1 − µ2 k22 dσ12
2
′ ′
+ 2 − d + 2d ln
dKL (p kq ) =
2
σ1
σ22
σ2
σ 
πkµ1 − µ2 k22 d(σ12 − σ22 )
2
+
+
d
ln
.
=
σ1
σ22
2σ22
Hence, by Pinsker inequality
′

′

kp − q k1 ≤
By triangle inequality,

p

2dKL

(p′ kq ′ )

≤

√

2πkµ1 − µ2 k2
+
σ2

p

d|σ12 − σ22 |
+
σ2

r

2d log

σ 
2
.
σ1

kp − qk1 ≤ kp − w2 p′ k1 + kw2 p′ − w2 q ′ k1 + kw2 q ′ − qk1 ≤ |w1 − w2 | + w2 kp′ − q ′ k1 ,
which gives the required bound. An identical proof works when w1 ≤ w2 .
47

Higher Dimensional Gaussians and Approximations. Let γd be the Gaussian measure
associated with a standard Gaussian with mean 0 and variance 1 in each direction.
Using concentration bounds for the χ2 random variables, we have the following bounds for the
lengths of vectors picked according to a standard Gaussian in d dimensions (see (4.3) in [31]).
Lemma A.4. For a standard Gaussian in d dimensions (mean 0 and variance 1/(2π) in each
direction), and any t > 0
i
h
√
1
(d + 2 dt + 2t) ≤ e−t .
P kxk2 ≥
x∼γd
2π
h
√ i
1
(d − 2 dt) ≤ e−t .
P kxk2 ≤
x∼γd
2π

Similarly, the following lemma shows a simple bound for the truncated moments, when x is
generated according to N (0, σ 2 /2π)d .
Lemma A.5. For any τ ≥ q, and any q ∈ Z+
Z
q
2
2
q
√ kxk2 exp(−πkxk2 /σ )dx ≤ σ exp(−4d).

(88)

kxk2 ≥2q dσ

√ √
Proof. Assume w.l.o.g that σ = 1. For kxk2 ≥ 2q d/ 2π, kxkq2 ≤ exp(πkxk22 /2). Hence,
Z
Z
q
2
2
kxk
exp(−πkxk
)dx
≤
2
√
√ exp(−πkxk /2)dx
2
kxk2 ≥2q d
kxk2 ≥2q d
Z


2
d/2
≤2
√ exp −πkyk2 /2 d y ≤ exp(−5d + d/2) ≤
kyk2 ≥ 2q√

d
2

1
,
16πd2

where y is distributed as a normal d-dimensional r.v. with mean 0 and variance 1 in each direction.
Fact A.6 (Stirling Approximation). For any n ≥ 1,

B

√

√
2πn(n/e)n ≤ n! ≤ e n(n/e)n .

Newton’s method for solving non-linear equations

We use a standard theorem that shows quadratic convergence of the Newton method in any normed
space [5] in the restricted setting where both the range and domain of f is Rm .
Consider a system of m non-linear equations in variables u1 , u2 , . . . , um :
∀j ∈ [m], fj (u1 , . . . , um ) = bj .
Newton’s method starts with an initial point u(0) close to a solution u∗ of the non-linear system.
Formally, u(0) ∈ N where N is an appropriately defined neighborhood set N = {y : ky − u∗ k ≤ ε0 }.
Let F ′ (u) = Jf (u) ∈ Rm×m be the Jacobian of the system given by the non-linear functional
∂fj (u)
|y is evaluated at
f : Rm → Rm , where the (j, i)th entry is the partial derivate Jf (j, i) = ∂u
i
y. Additionally, for our algorithm, we assume that given any y ∈ N , we only have access to an
estimates b̃, F̃ (u), F̃ ′ (u) of vector b ∈ Rm , F (u) ∈ Rm and F ′ (u) ∈ Rm×m respectively 8 .
8

These errors in the estimate may occur due to sampling errors or precision errors.

48

Newton’s method starts with the initializer u(0) , and updates the solution using the iteration:

−1 

u(t+1) = u(t) + F̃ ′ (u(t) )
b̃ − f (u(t) ) .
(89)

The convergence error will be measured in the ℓp norm for any p ≥ 1. In what follows, kxk :=
kxkp for x ∈ Rm , kM k := kM kp→p for M ∈ Rm×m and kT k := kT kp×p→p for T ∈ Rm×m×m . We
first state a simple mean-value theorem that will be useful in the analysis of the Newton method,
as well as in applying the guarantees in the context of mixtures of Gaussians.
Lemma B.1 (Proposition 5.3.11 in [5]). Consider a function H : K ⊂ Rm1 → Rm2 , with K being
an open set. Assume H is differentiable on K and that F ′ (u) is a continuous function of u on K.
Assume u, w ∈ K and assume line segment from joining them is also contained in K. Then
kF (u) − F (w)k ≤ sup kF ′ ((1 − θ)u + θw)k · ku − wk.
0≤θ≤1

The following theorem gives robustness guarantees for Newton’s method. It is obtained by using
matrix perturbation analysis along with a standard theorem regarding the quadratic convergence
of the Newton’s method (see Theorem 5.4.1 in [5]).
Theorem B.2. Assume u∗ ∈ Rm is a solution to the equation F (u) = b where F : Rm → Rm such
that J −1 = (F ′ )−1 exists in a neighborhood N = {y : ky − u∗ k ≤ ku(0) − u∗ k}, and F ′ : Rm → Rm×m
is locally L-Lipschitz continuous in the neighborhood N i.e.,
kF ′ (u) − F ′ (v)k ≤ Lku − vk ∀u, v ∈ N .
Further, let the estimates b̃, F̃ (u), F̃ ′ (u) satisfy for some η1 , η2 , η3 > 0 and all y ∈ N
kb̃ − bk ≤ η1 , kF̃ (u) − F (u)k ≤ η2 , kF̃ ′ (u) − F ′ (u)k ≤ η3 .
Then if η3 kF ′ (u(t) )−1 k < 1, kF ′ (u)k ≤ B, then for all u ∈ N , the error εt = ku(t) − u∗ k after the t
iterations of (89) satisfies


εt+1 ≤ ε2t · LkF ′ (u(t) )−1 k + kF ′ (u(t) )−1 k η1 + η2 + 4η3 εt kF ′ (u(t) )−1 kB
(90)

Proof. We sketch the proof here. For convenience, let us use A := F ′ (u(t) ) and Ã := F̃ ′ (u(t) ) to
denote the derivate and its estimate at u(t) . Let z (t+1) be the Newton iterate after the (t + 1)th
step if we had the exact values of F (u(t) and F ′ (u(t) ) i.e., z (t+1) = u(t) + A−1 b − F (u(t) ) . From
the standard analysis of the Newton method (see Theorem 5.4.1 in [5]),
kz (t+1) − u∗ k ≤ LkA−1 kku(t) − u∗ k2 .

Further, the error between the actual Newton update u(t+1) and z (t) due to the estimates Ã and
F̃ (u(t) ) is




u(t+1) − z (t+1) = Ã−1 b̃ − F̃ (u(t) ) − A−1 b − F (u(t) )




= (Ã−1 − A−1 ) b̃ − F̃ (u(t) ) + A−1 b̃ − b + F̃ (u(t) ) − F (u(t) )
ku(t+1) − z (t+1) k ≤ kÃ−1 − A−1 kkb̃ − F̃ (u(t) )k + kA−1 k(η1 + η2 ).
49

From perturbation bounds on matrix inverses [13], if kA−1 Ek < 1,
kA−1 − (A + E)−1 k ≤ kA−1 k ·

kA−1 Ek
.
1 − kA−1 Ek

Also from Lemma B.1, kb−F (u(t) )k ≤ kF ′ (u′ )kku(t) −u∗ k for some u′ ∈ N . Substituting E = Ã−A,

kA−1 (Ã − A)k 
η1 + η2 + kb − F (u(t) )k + kA−1 k(η1 + η2 )
1 − kA−1 (Ã − A)k


≤ 2η3 kA−1 k2 η1 + η2 + kF ′ (u′ )kku(t) − u∗ k + kA−1 k(η1 + η2 )

ku(t+1) − z (t+1) k ≤ kA−1 k

≤ 4η3 εt kA−1 k2 kF ′ (u′ )k + kA−1 k(η1 + η2 ).

Remark B.3. While the above theorem requires that the derivative F ′ is locally L-Lipschitz, this
is a weaker condition than requiring a upper bound on the operator norm of the second derivative
F ′′ . Lemma B.1 shows that it also suffices if kF ′′ (u)k ≤ L for all u ∈ N .
Corollary B.4. Under the conditions of Theorem B.2, there exists 0 < ε0 <
that for any given δ ∈ (0, 1), there is an η1 , η2 , η3 > 0 with
(η1 + η2 ) <

δ
4kF ′ (u(t) )−1 k

and η3 <

δ
4kF ′ (u(t) )−1 k2

· min{1,

1
,
2LkF ′ (u(t) )−1 k

such

1
}
B

such that after T = log log(1/δ) iterations of the Newton’s method, we have
ku(T ) − u∗ k ≤ δ.
Proof. For the given setting of η1 , η2 , η3 , we have (η1 +η2 )kF ′ (u(t) )−1 k < 4δ , and BkF ′ (u(t) )−1 k2 εt η3 ≤
δ/4. From Theorem B.2, we have that for any t,
δ
εt+1 ≤ ε2t LkF ′ (u(t) )−1 k + .
2
Further ε1 ≤ ε20 LkF ′ (u(t) )−1 k <

1
.
4LkF ′ (u(t) )−1 k

By induction, it follows that
t

εt ≤

2−2
.
LkF ′ (u(t) )−1 k

Hence, this gives the required guarantee.

C

Dimension Reduction using PCA

Here we give a proof of the assertion that for mixtures of spherical Gaussians, we can assume
without loss of generality that d ≤ k.

50

Theorem C.1 (Same as Theorem 4.2). Let {(wi , µi , σi ) : i ∈ [k]} be a mixture of k spherical
Gaussians that is ρ bounded, and let wmin be the smallest mixing weight. Let µ′1 , µ′2 , . . . , µ′k be the
projections onto the subspace spanned by the top k singular vectors of sample matrix X ∈ Rd×N .
−1
For any ε > 0, with N = poly(d, ρ, wmin
, ε−1 ) samples we have with high probability
∀i ∈ [k], kµi − µ′i k2 ≤ ε.
Proof. Let δ = wmin ε4 /(2ρ2 ) and η = δ2 /2. Let A be the population average, i.e.,
X
1 X
wi σi2 .
A = E[xxT ] = M + σ̄ 2 I, where M =
wi µi µTi , and σ̄ 2 =
2π
i

i∈[k]

Let λ1 ≥ λ2 ≥ · · · ≥ λd ≥ 0 be the eigenvalues of M sorted in non-increasing order. Since M is
of rank at most k, we have λk+1 = · · · = λd = 0. Let r ≤ k be defined as the smallest index such
that λr+1 < δ. Let U be represent the orthogonal projector onto the top-r eigenspace of M (and
hence of A too), and U ⊥ = I − U . Notice that kU ⊥ M U ⊥ k = λr+1 < δ. Then, using the positive
semidefinite inequality wi U ⊥ µi µTi U ⊥  U ⊥ M U ⊥ , we obtain that
kU ⊥ µi k22 ≤ wi−1 U ⊥ M U ⊥ = λr+1 /wi < δ/wmin .
Next, let Ã represent the sample average Ã = N1 XX T (so Ã converges to A as N → ∞). Since
N ≥ poly(d)/η 2 , using standard concentration bounds (see Theorem 6.1.1 in [43] and related notes)
we have with high probability that kA − Ãk < η. Let Ũ be the orthogonal projector onto the top-k
eigenspace of Ã, and let Ṽ = I − Ũ be the orthogonal projector onto the bottom (d − k) eigenspace
of Ã. Hence for each i ∈ [k], µ′i = Ũ µi . Notice that from Weyl’s perturbation bounds for eigenvalues
(see [13], Theorem III.2.1) λk+1 (Ã) ≤ σ̄ 2 + η. Therefore, the eigenvalues of Ã corresponding to Ṽ
(which are all at most σ̄ 2 + η), and the eigenvalues of A corresponding to U (which are all at least
σ̄ 2 + δ) are separated by at least δ − η. From standard perturbation bounds for eigenvectors (see
[13], Theorem VII.3.1), we have
kU Ṽ k ≤

η
kA − Ãk
≤
≤ δ.
δ−η
δ−η

Hence, for each i ∈ [k],
kµi − µ′i k22 = kṼ µi k22 = hµi , Ṽ µi i

= hU µi , Ṽ µi i + hU ⊥ µi , Ṽ µi i = hµi , U Ṽ µi i + hU ⊥ µi , Ṽ µi i

≤ kU Ṽ kkµi k22 + kU ⊥ µi k2 kµi k2
r
δ
2
≤ δkµi k2 +
· kµi k2 ≤ ε2
wmin
by our choice of δ = wmin ε4 /(2ρ2 ).

Acknowledgements
The authors thank Santosh Vempala for suggesting the problem of learning one-dimensional Gaussians, and other helpful discussions. Part of this work was done when the second author was at the
Courant Institute and the Simons Collaboration on Algorithms and Geometry.
51

References
[1] Emmanuel Abbe. Community detection and the stochastic block model. 2016. Available from
http://www.ee.princeton.edu/research/eabbe.
[2] Dimitris Achlioptas and Frank McSherry. On spectral learning of mixtures of distributions.
In Learning Theory, pages 458–469. Springer, 2005.
[3] Joseph Anderson, Mikhail Belkin, Navin Goyal, Luis Rademacher, and James R. Voss. The
more, the merrier: the blessing of dimensionality for learning large Gaussian mixtures. In
Proceedings of The 27th Conference on Learning Theory, COLT 2014, Barcelona, Spain, June
13-15, 2014, pages 1135–1164, 2014.
[4] Sanjeev Arora and Ravi Kannan. Learning mixtures of arbitrary Gaussians. In Proceedings of
the thirty-third annual ACM symposium on Theory of computing, pages 247–257. ACM, 2001.
[5] Kendall Atkinson and Weimin Han. Theoretical numerical analysis : a functional analysis
framework. Texts in applied mathematics. Springer, New York, Berlin, Paris, etc., 2001.
[6] Pranjal Awasthi, Avrim Blum, and Or Sheffet. Center-based clustering under perturbation
stability. Information Processing Letters, 112(12):49 – 54, 2012.
[7] Pranjal Awasthi and Or Sheffet. Improved spectral-norm bounds for clustering. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages
37–49. Springer, 2012.
[8] Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu. Statistical guarantees for the
EM algorithm: From population to sample-based analysis. CoRR, abs/1408.2156, 2014.
[9] Mikhail Belkin and Kaushik Sinha. Polynomial learning of distribution families. In Foundations
of Computer Science (FOCS), 2010 51st Annual IEEE Symposium on, pages 103–112. IEEE,
2010.
[10] Aditya Bhaskara, Moses Charikar, Ankur Moitra, and Aravindan Vijayaraghavan. Smoothed
analysis of tensor decompositions. In Proceedings of the 46th Symposium on Theory of Computing (STOC). ACM, 2014.
[11] Aditya Bhaskara, Moses Charikar, and Aravindan Vijayaraghavan. Uniqueness of tensor decompositions with applications to polynomial identifiability. Proceedings of the Conference on
Learning Theory (COLT)., 2014.
[12] Aditya Bhaskara, Ananda Suresh, and Morteza Zadimoghaddam. Sparse Solutions to Nonnegative Linear Systems and Applications. In Guy Lebanon and S. V. N. Vishwanathan, editors,
Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics,
volume 38 of Proceedings of Machine Learning Research, pages 83–92, San Diego, California,
USA, 09–12 May 2015. PMLR.
[13] Rajendra Bhatia. Matrix Analysis, volume 169. Springer, 1997.

52

[14] Spencer Charles Brubaker and Santosh Vempala. Isotropic pca and affine-invariant clustering.
In Proceedings of the 2008 49th Annual IEEE Symposium on Foundations of Computer Science,
FOCS ’08, pages 551–560, Washington, DC, USA, 2008. IEEE Computer Society.
[15] Siu-On Chan, Ilias Diakonikolas, Rocco A. Servedio, and Xiaorui Sun. Efficient density estimation via piecewise polynomial approximation. In Proceedings of the Forty-sixth Annual
ACM Symposium on Theory of Computing, STOC ’14, pages 604–613, New York, NY, USA,
2014. ACM.
[16] Sanjoy Dasgupta. Learning mixtures of Gaussians. In Foundations of Computer Science, 1999.
40th Annual Symposium on, pages 634–644. IEEE, 1999.
[17] Sanjoy Dasgupta and Leonard Schulman. A probabilistic analysis of EM for mixtures of
separated, spherical Gaussians. The Journal of Machine Learning Research, 8:203–226, 2007.
[18] Constantinos Daskalakis and Gautam Kamath. Faster and sample near-optimal algorithms
for proper learning mixtures of Gaussians. In Maria-Florina Balcan, Vitaly Feldman, and
Csaba Szepesvári, editors, Proceedings of The 27th Conference on Learning Theory, COLT
2014, Barcelona, Spain, June 13-15, 2014, volume 35 of JMLR Workshop and Conference
Proceedings, pages 1183–1213. JMLR.org, 2014.
[19] Constantinos Daskalakis, Christos Tzamos, and Manolis Zampetakis. Ten steps of EM suffice
for mixtures of two Gaussians. CoRR, abs/1609.00368, 2016.
[20] Ilias Diakonikolas, Daniel M. Kane, and Alistair Stewart. Statistical query lower bounds for robust estimation of high-dimensional Gaussians and Gaussian mixtures. CoRR, abs/1611.03473,
2016.
[21] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.13 of 201609-16. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert,
C. W. Clark, B. R. Miller and B. V. Saunders, eds.
[22] John
Duchi.
Derivations
for
linear
algebra
and
web.stanford.edu/~jduchi/projects/general_notes.pdf. [Online].

optimization.

[23] Jon Feldman, Rocco A. Servedio, and Ryan O’Donnell. PAC learning axis-aligned mixtures
of Gaussians with no separation assumption. In Proceedings of the 19th annual conference on
Learning Theory, COLT’06, pages 20–34, Berlin, Heidelberg, 2006. Springer-Verlag.
[24] 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, STOC 2015, Portland, OR, USA, June 14-17, 2015, pages 761–770, 2015.
[25] Navin Goyal, Santosh Vempala, and Ying Xiao. Fourier PCA and robust tensor decomposition.
In Symposium on Theory of Computing, STOC 2014, New York, NY, USA, May 31 - June
03, 2014, pages 584–593, 2014.
[26] Moritz Hardt and Eric Price. Tight bounds for learning a mixture of two Gaussians. In
Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC
2015, Portland, OR, USA, June 14-17, 2015, pages 753–760, 2015.
53

[27] 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.
[28] Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant. Efficiently learning mixtures of
two Gaussians. In Proceedings of the 42nd ACM symposium on Theory of computing, pages
553–562. ACM, 2010.
[29] Ravindran Kannan, Hadi Salmasian, and Santosh Vempala. The spectral method for general
mixture models. SIAM J. Comput., 38(3):1141–1156, 2008.
[30] Amit Kumar and Ravindran Kannan. Clustering with spectral norm and the k-means algorithm. In Foundations of Computer Science (FOCS), 2010 51st Annual IEEE Symposium on,
pages 299–308. IEEE, 2010.
[31] B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection.
Ann. Statist., 28(5):1302–1338, 10 2000.
[32] Laurent Massoulié. Community detection thresholds and the weak ramanujan property. In
Proceedings of the 46th Annual ACM Symposium on Theory of Computing, STOC ’14, pages
694–703, New York, NY, USA, 2014. ACM.
[33] F. McSherry. Spectral partitioning of random graphs. In Proceedings of the 42nd IEEE
symposium on Foundations of Computer Science, FOCS ’01, pages 529–, Washington, DC,
USA, 2001. IEEE Computer Society.
[34] D. G. Mixon, S. Villar, and R. Ward. Clustering subgaussian mixtures with k-means. In 2016
IEEE Information Theory Workshop (ITW), pages 211–215, Sept 2016.
[35] Ankur Moitra and Gregory Valiant. Settling the polynomial learnability of mixtures of Gaussians. In Foundations of Computer Science (FOCS), 2010 51st Annual IEEE Symposium on,
pages 93–102. IEEE, 2010.
[36] Elchanan Mossel, Joe Neeman, and Allan Sly. A proof of the block model threshold conjecture.
CoRR, abs/1311.4115, 2013.
[37] Elchanan Mossel, Joe Neeman, and Allan Sly. Belief propagation, robust reconstruction and
optimal recovery of block models. In Proceedings of the Conference on Learning Theory, pages
356–370, 2014.
[38] Karl Pearson. Contributions to the mathematical theory of evolution. Philosophical Transactions of the Royal Society of London. A, 185:71–110, 1894.
[39] Nathan Srebro, Gregory Shakhnarovich, and Sam Roweis. An investigation of computational
and informational limits in Gaussian mixture clustering. In Proceedings of the 23rd International Conference on Machine Learning, ICML ’06, pages 865–872, New York, NY, USA, 2006.
ACM.
[40] Ananda Theertha Suresh, Alon Orlitsky, Jayadev Acharya, and Ashkan Jafarpour. Nearoptimal-sample estimators for spherical Gaussian mixtures. In Z. Ghahramani, M. Welling,
54

C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, Advances in Neural Information
Processing Systems 27, pages 1395–1403. Curran Associates, Inc., 2014.
[41] Henry Teicher. Identifiability of mixtures. The annals of Mathematical statistics, 32(1):244–
248, 1961.
[42] Henry Teicher. Identifiability of mixtures of product measures. The Annals of Mathematical
Statistics, 38(4):1300–1302, 1967.
[43] Joel A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathematics, 12(4):389–434, Aug 2012.
[44] J.M. Varah. A lower bound for the smallest singular value of a matrix. Linear Algebra and its
Applications, 11(1):3 – 5, 1975.
[45] Santosh Vempala and Grant Wang. A spectral algorithm for learning mixture models. Journal
of Computer and System Sciences, 68(4):841–860, 2004.
[46] Ji Xu, Daniel J. Hsu, and Arian Maleki. Global analysis of expectation maximization for
mixtures of two Gaussians. In NIPS, 2016.

55

