Settling the Polynomial Learnability of Mixtures of Gaussians

Introduction

Given access to random samples generated from a mixture of (multivariate) Gaussians, the algorithmic problem of learning the parameters of the underlying distribution is of fundamental importance in physics, biology, geology, social sciences – any area in which such finite mixture models arise . Starting with Dasgupta , a series of work in theoretical computer science has sought to find (or disprove the existence of) an efficient algorithm for this task . In this paper, we settle the polynomial-time learnability of mixtures of Gaussians, giving an algorithm that uses a polynomial amount of data and estimates the components at an inverse polynomial rate under provably minimal assumptions on the mixture (specifically, that the mixing weights and the statistical distance between the components are bounded away from zero). As a corollary, our efficient learning algorithm can be employed to yield the first provably efficient algorithm for near-optimal clustering and density estimation, without any restrictions on the Gaussian mixture. Finally, we note that the runtime and data requirements of our algorithm are exponential in the number of Gaussian components; however, as we show in Section 6, this exponential dependence is necessary. In the remainder of this section, we briefly summarize previous work on this problem, formally state our main result, and then discuss the differences between learning mixtures of $2$ Gaussians, and mixtures of many Gaussians, which motivates the high-level outline of our algorithm presented in Section 2. We first define a Gaussian Mixture Model (GMM).

The most popular solution for recovering reasonable estimates of the components of GMMs in practice is the EM algorithm given by Dempster, Laird and Rubin . This algorithm is a local-search heuristic that converges to a set of parameters that locally maximizes the probability of generated the observed samples. However, the EM algorithm is a heuristic only, and makes no guarantees about converging to an estimate that is close to the true parameters. Worse still, the EM algorithm (even for univariate mixtures of just two Gaussians) has been observed to converge very slowly (see Redner and Walker for a thorough treatment ).

In order to even hope for an algorithm (not necessarily even polynomial time), we would need a uniqueness property – that two distinct mixtures of Gaussians must have different probability density functions. Teicher demonstrated that a mixture of Gaussians can be uniquely identified (up to a relabeling components) by considering the probability density function at points sufficiently far from the centers (in the tails). However, such a result sheds little light on the rate of convergence of an estimator: If distinguishing Gaussian mixtures really required analyzing the tails of the distribution, then we would require an enormous number of data samples!

Dasgupta introduced theoretical computer science to the algorithmic problem of provably recovering good estimates for the parameters in polynomial time (and a polynomial number of samples). His technique is based on projecting data down to a randomly chosen low-dimensional subspace, finding an accurate clustering. Given enough accurately clustered points, the empirical means and co-variances of these points will be a good estimate for the actual parameters. Arora and Kannan extended these ideas to work in the much more general setting in which the co-variances of each Gaussian component could be arbitrary, and not necessarily almost spherical as in . Yet both of these techniques are based on the concentration of distances (under random projections), and consequently required that the centers of the components be separated by at least $\sqrt{n}$ times the largest variance. Vempala and Wang and Achlioptas and McSherry introduced the use of spectral techniques, and were able to overcome this barrier (of relying on distance concentration) by choosing a subspace on which to project based on large principle components. Brubaker and Vempala later gave the first affine-invariant algorithm for learning mixtures of Gaussians, and these ideas proved to be central in subsequent work .

Yet all of these approaches for provably learning good estimates require, at the very least, that the statistical overlap (i.e. one minus the statistical distance) between each pair of components be at least smaller than some constant (in some cases, it is even required that the statistical overlap be exponentially small). Recently, Felman et al gave a polynomial time algorithm for the related problem of density estimation (without any separation condition) for the special case of axis-aligned GMMs (GMMs where each component has principle coordinates aligned with the coordinate axes). Also without any separation requirements, Belkin and Sinha showed that one can efficiently learn GMMs in the special case that all components are identical spherical Gaussians. Most similar to the present work is the recent work of Kalai et al , that gave a learning algorithm for the case of mixtures of two arbitrary Gaussians with provably minimal assumptions.

2 Main Results

In this section we state our main results. To motivate these results, we first state three obvious lower bounds for recovering the parameters of a GMM $F=\sum_{i=1}^{k}w_{i}F_{i},$ which motivate our defintion of $\epsilon$-statistically learnable. We provide a formal definition of statistical distance in Section 2.1.

Permuting the order of the components does not change the resulting density, thus at best the hope is to recover the parameter set, $\{(w_{1},\mu_{1},\Sigma_{1}),\ldots,(w_{k},\mu_{k},\Sigma_{k})\}.$

We require at least $\Omega(1/\min_{i}(w_{i}))$ samples to estimate the parameters, since we require this number of samples to ensure that we have seen, with reasonable probability, any sample from each component.

If $F_{i}=F_{j}$, then it is impossible to accurately estimate $w_{i},$ and in general we require at least $\Omega(1/D(F_{i},F_{j}))$ samples to estimate $w_{i}$, where $D(F_{i},F_{j})$ denotes the statistical distance between the two distributions.

We call a GMM $F=\sum_{i}w_{i}F_{i}$ $\epsilon$-statistically learnable if $\min_{i}w_{i}\geq\epsilon$ and $\min_{i\neq j}D(F_{i},F_{j})\geq\epsilon$.

We now consider what it means to “accurately recover the mixture components”.

Given two $n$-dimensional GMMs of $k$ Gaussians, $F=\sum_{i}w_{i}\mathcal{N}(\mu_{i},\Sigma_{i})$ and $\hat{F}=\sum_{i}\hat{w}_{i}\mathcal{N}(\hat{\mu}_{i},\hat{\Sigma}_{i})$, we call $\hat{F}$ an $\epsilon$-close estimate for $F$ if there is permutation function $\pi:[k]\rightarrow[k]$ such that for all $i\in[k]$

$D(\mathcal{N}(\mu_{i},\Sigma_{i}),\mathcal{N}(\hat{\mu}_{\pi(i)},\hat{\Sigma}_{\pi(i)}))\leq\epsilon$

Note that the above definition of an $\epsilon$-close estimate is affine invariant. This is more natural than defining a good estimate in terms of additive errors, since in general, even estimating the mean of an arbitrary Gaussian to some fixed additive precision is impossible without restrictions on the covariance, as scaling the data will scale the error linearly. We can now state our main theorem:

Given any $n$ dimensional mixture of $k$ Gaussians $F$ that is $\epsilon$-statistically learnable, we can output an $\epsilon$-close estimate $\hat{F}$ and the running time and data requirements of our algorithm (for any fixed $k$) are polynomial in $n$, and $\frac{1}{\epsilon}$.

The guarantee in the main theorem implies that the estimated parameters are off by an additive $O(\epsilon\sigma^{2}_{max})$, where $\sigma^{2}_{max}$ is the largest (projected) variance of any Gaussian in any direction.

Throughout this paper, we favor clarity of proof and exposition above optimization of runtime. Since our main goal is show that these problems can be solved in polynomial time, we make very little effort to optimize the exponent. Our algorithms are polynomial in the dimension, inverse of the success probability, and inverse of the target accuracy for any fixed number of Gaussians, $k$. The dependency on $k$, however, is severe: the degree of our polynomials are linear in $k$. In Section 6, we give a natural construction of two GMMs $F,F^{\prime}$ of $k$ components that are each $1/k$-statistically learnable, satisfy $D(F,F^{\prime})\leq e^{-k}$, but $F$ is not even a $1/4$-close estimate of $F$. Thus we require an exponential in $k$ number of samples to even distinguish these two mixtures, demonstrating that the exponential dependency on $k$ in our learning algorithms is inevitable.

There exists two GMMs $F,F^{\prime}$ of $k$ components each that satisfies the following properties:

$F,F^{\prime}$ are $1/k$-statistically learnable.

$F$ is not a $1/4$-close estimate of $F^{\prime}$.

3 Applications

We can leverage our main theorem to show that we can efficiently perform density estimation for arbitrary GMMs. For density estimation—as opposed to parameter recovery—we only care to recover a distribution that is similar to the GMM, without worrying about matching each component; in particular, if the true weight of one of the components is negligible, we can simply disregard that component with negligible effect on the statistical distance; if two components are nearly identical in statistical distance, we can simply regard them as being merged into one component. For these reasons, we can perform density estimation efficiently without the restriction to $\epsilon$-statistically learnable distributions, that was required for Theorem 1.

For any $n\geq 1,$ $\epsilon,\delta>0,$ and any $n$-dimensional GMM $F=\sum_{i=1}^{k}w_{i}F_{i},$ given access to independent samples from $F$, there is an algorithm that outputs $\hat{F}=\sum_{i=1}^{k}\hat{w_{i}}\hat{F_{i}}$ such that with probability at least $1-\delta$ over the randomization in the algorithm and in selecting the samples, $D(F,\hat{F})\leq\epsilon.$ Additionally, the runtime and number of samples is bounded by $poly(n,1/\epsilon,1/\delta).$

The proof of this corollary follows immediately from combining our main theorem, with the arguments in Appendix D. In fact, an almost identical approach to how we construct the General Univariate Algorithm from the Basic Univariate Algorithm (again in Appendix D) will work because we can run our main algorithm with many different parameter ranges so that most estimates are correct, and determine a consensus among the estimate so that we can recover a good statistical approximation to $F$ without any assumptions on the mixture - not even $\epsilon$-statistical learnability.

For any $n\geq 1,$ $\epsilon,\delta>0,$ and any $n$-dimensional $\epsilon$-statistically learnable GMM $F=\sum_{i=1}^{k}w_{i}F_{i},$ given access to independent samples from $F$, there is an algorithm that outputs a classifier $C^{\prime}_{F}$ such that with probability at least $1-\delta$ over the randomization in the algorithm and in selecting the samples, the error of $C_{F}$ is at most $\epsilon$ larger than the error of any classifier $C^{\prime}$. Additionally, the runtime and number of samples used is bounded by $poly(n,1/\epsilon,1/\delta).$

The proof of this corollary follows immediately from our main theorem (yet here we need the assumption of $\epsilon$-statistical learnability in this case).

4 Comparing Learning Two Gaussians to Learning Many

This work leverages several key ideas initially presented in which were used to show that learning mixtures of two arbitrary Gaussians can be done efficiently. Nevertheless, additional high-level insights, and technical details were required to extend the previous work to give an efficient learning algorithm for an arbitrary mixture of many Gaussians. In this section we briefly summarize the algorithm for learning mixtures of two Gaussians given in , and then describe the hurdles to extending it to the general case. This discussion will provide insights and motivate the high-level structure of the algorithm presented in this paper, as well as clarify which components of the proof are new, and which are straight-forward adaptations of ideas from .

Throughout this discussion, it will be helpful to refer to parameters $\epsilon_{1},\epsilon_{2},\epsilon_{3},$ which are polynomially related to each other, and satisfy $\epsilon_{1}<<\epsilon_{2}<<\epsilon_{3}$.

There are three key components to the proof that mixtures of two Gaussians can be learned efficiently: the 1-d Learnability Lemma, the Random Projection Lemma, and the Parameter Recovery Lemma. The 1-d Learnability Lemma states that given a mixture of two univariate Gaussians whose two components have nonnegligible statistical distance, one can efficiently recover accurate estimates of the parameters of the mixture. It is worth noting that in the univariate case, saying that the statistical distance between two Gaussians is non-negligible is roughly equivalent (polynomially related) to saying that the two sets of parameters are non-negligibly different, ie. the parameter distance, $|\mu-\mu^{\prime}|+|\sigma^{2}-\sigma^{\prime 2}|,$ is non-negligible. The Random Projection Lemma states that, given an $n$-dimensional mixture of two Gaussians which is in isotropic position and whose components have nonnegligible statistical distance, with high probability over the choice of a random unit vector $r,$ the projection of the mixture onto $r$ will yield a univariate mixture of two Gaussians that have nonnegligible statistical distance (say $\epsilon_{3}$). The final component—the Parameter Recovery Lemma—states that, given a Gaussian $G$ in $n$ dimensions, if one has extremely accurate estimates (say to within some $\epsilon_{1}$) of the mean and variance of $G$ projected onto $n^{2}$ sufficiently distinct directions (directions that differ by at least $\epsilon_{2}>>\epsilon_{1}$) one can accurately recover the parameters of $G$.

Given these three pieces, the high-level algorithm for learning mixtures of two Gaussians is straight-forward:

Pick $n^{2}$ vectors $r_{1},\ldots,r_{n^{2}},$ that are “close” to $r$, say $|r_{i}-r|\approx\epsilon_{2}.$

For each $i=1,\ldots,n^{2},$ learn extremely accurate (to accuracy $\epsilon_{1}<<\epsilon_{2}$) univariate parameters $w_{i},\mu_{i},\sigma_{i},\mu^{\prime}_{i},\sigma^{\prime}_{i}$ for the projection of the mixture onto the vector $r_{i}$.

Since $|r_{i}-r_{j}|\approx\epsilon_{2},$ it is not hard to show that with high probability, $|\mu_{i}-\mu_{j}|<<\epsilon_{3},|\sigma_{i}-\sigma_{j}|<<\epsilon_{3}$ and by the Random Projection Lemma, $||(\mu_{i},\sigma_{i})-(\mu^{\prime}_{i},\sigma^{\prime}_{i})||>>\epsilon_{3}$ thus it will be easy to accurately match up which parameters come from which component in the different projections, and we can apply the Parameter Recovery Lemma to each of the two components.

Some of the above ideas are immediately applicable to the problem of learning mixtures of many Gaussians: we can clearly use the Parameter Recovery Lemma without modification. Additionally, we prove a generalization of the 1-d Learnability Lemma for mixtures of arbitrary numbers of Gaussians, provided each component has non-negligible statistical distance (which, while technically tedious, employs the key idea from of “deconvolving” by a suitably chosen Gaussian—see Appendix B). Given this extension, if we were given a mixture of $k$ Gaussians in isotropic position, and were guaranteed that the projection onto some vector $r$ resulted in a univariate mixture of Gaussians for which all pairs of components either had reasonably different means or reasonably different variances, then we could piece together the parts more-or-less as in the 2-Gaussians case.

Unfortunately, however, the Random Projection Lemma, ceases to hold in the general setting. There exist mixtures of just three Gaussians with significant pairwise statistical distances, that are in isotropic position, but have the property that with extremely high probability over choices of random unit vector $r$, the projection of the mixture onto $r$ yields a distribution that is extremely close to a univariate mixture of two Gaussians. This observation would foil the approach employed in the case of just two Gaussians! Another difficulty is that if we take $n^{2}$ slightly different projections of our mixture of $k$ Gaussians, then it is possible that in some of the projections we see what looks like a mixture of $k^{\prime}<k$ univariate Gaussians, and in some other projections we see what looks like a mixture of $k^{\prime\prime}$ univariate Gaussians. How do we match up estimates from projections onto different directions when the number of Gaussians in the estimate can differ? Or what if each projection results in an estimate that is a mixture of $k^{\prime}<k$ Gaussians. Then how can we recover an $n$-dimensional estimate that is a mixture of $k$ Gaussians?

Outline and Definitions

We now discuss the high-level structure of our learning algorithm, building from the intuition given in the preceding section. At the highest level, our learning algorithm has the following form: Given access to samples from a mixture of $k$ Gaussians,

Learn the parameters of some mixture of $k^{\prime}\leq k$ Gaussians, where each learned Gaussian component roughly corresponds to one or more of the Gaussians in the original mixture.

If $k^{\prime}<k$, for each of the $k^{\prime}$ components recovered in the previous step, examine it closely and figure out whether it corresponds to a single Gaussian component of the original mixture, or whether it is a mixture of several of the original components (in which case we will then need to learn the parameters of these sub-components).

To accomplish the first step, we will require accurate parameters of the projection of each of the $k^{\prime}$ “clusters” of components, onto $n^{2}$ univariate projections. To do this, we employ a robust univariate algorithm which, given access to samples from a univariate GMM, essentially searches for some target resolution window $(w_{1},w_{2})$ with $w_{1}<<w_{2},$ such that the GMM is very close ($w_{1}$-close) to a GMM of $k^{\prime}\leq k$ statistically very distinct components (each pair of components is at least $w_{2}$ far apart).

Given our robust univariate algorithm, we embark on a partition pursuit where we try to find $n^{2}$ vectors that yield consistent and compatible univariate parameter sets–in particular, we require that each of the $n^{2}$ univariate projections yields parameters that satisfy three conditions: 1) they have the same number of components, 2) the recovered parameters are much more precise than the distances between the $n^{2}$ projections, and 3) that the distance between the components is large enough so as to ensure an accurate matching of the components in the different projections.

Finally, given the ability to accurately recover $k^{\prime}\leq k$ high-dimensional Gaussians, where each learned Gaussian component roughly corresponds to one or more of the Gaussians in the original mixture, we want to be able to examine each recovered component, and determine whether it corresponds to a single component of the original mixture, or a set of original components. We first claim that, with high probability, the only way a subset of original components will end up being grouped into a single recovered component is if the covariance of the mixture of that subset of components has a very small minimum eigenvalue. The existence of such an eigenvalue implies that we can accurately cluster the given sample points (whose covariance, recall, is roughly 1). Thus, given a recovered set of $k^{\prime}<k$ parameters, we examine one of these $k^{\prime}$ components; if the minimum eigenvalue is sufficiently small, we project the set of data samples onto the corresponding eigenvector, and then partition the sample points into two clusters (provided the eigenvalue is sufficiently small, since the overall mixture is in roughly isotropic position, we cluster so as to almost exactly respect some partition of the original components). Given the set of sample points corresponding (roughly) to the recovered component that had small eigenvalue, we simply re-scale the data so that this subsample is now in isotropic position, and recursively run the entire algorithm on this rescaled subsample of the data, which, as we argue, consists of a mixture of $k^{\prime\prime}<k$ components of the original mixture, with high probability. We call this clustering step hierarchical clustering.

We give a detailed summary in Appendix A of the main elements of each of these three main components: the robust univariate algorithm, partition pursuit, and hierarchical clustering.

Given two probability distributions $f(x),g(x)$ on $\Re^{n}$ we can define the statistical distance between these distributions as

We will also be interested in a related notion of the parameter distance between two univariate Gaussians:

Given two univariate Gaussians, $F_{1}=\mathcal{N}(\mu_{1},\sigma_{1}^{2}),F_{2}=\mathcal{N}(\mu_{2},\sigma_{2}^{2})$ we define the parameter distance as

In general, the parameter distance and the statistical distance between two univariate Gaussians can be unrelated. There are pairs of univariate Gaussians with arbitrarily small parameter distance, and yet statistical distance close to $1$, and there are pairs of univariate Gaussians with arbitrarily small statistical distance, and yet arbitrarily large parameter distances. But these scenarios can only occur if the variances can be arbitrarily small or arbitrarily large. In many instances in this paper, we will have reasonable upper and lower bounds on the variances and this will allow us to move back and forth from statistical distance and parameter distance, but we will highlight when we are doing so and note why we are able to assume an upper and lower bound on variance in that particular situation.

As we noted, there are $\epsilon$-statistically learnable mixtures of three Gaussians that are in isotropic position, but for which with overwhelming probability over a random direction $r$, in the projection onto $r$, there will be some pair of univariate Gaussians that are arbitrarily close in parameter distance. In these cases, our univariate algorithm may not return an estimate with three components, but will return a mixture which has only two components but is still a good estimate for the parameters of the projected mixture. To formalize this notion, we introduce what we call an $\epsilon$-correct sub-division.

Given a GMM of $k$ Gaussians, $F=\sum_{i}w_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2})$ and a GMM of $k^{\prime}\leq k$ Gaussians $\hat{F}=\sum_{i}\hat{w}_{i}\mathcal{N}(\hat{\mu}_{i},\hat{\sigma}_{i}^{2})$, we call $\hat{F}$ an $\epsilon$-correct subdivision of $F$ if there is a function $\pi:[k]\rightarrow[k^{\prime}]$ that is onto and

$\forall_{j\in[k^{\prime}]}|\sum_{i|\pi(i)=j}w_{i}-\hat{w}_{j}|\leq\epsilon$

$\forall_{i\in[k]}D_{p}(F_{i},\hat{F}_{\pi(i)})\leq\epsilon$

When considering high-dimensional mixtures, we replace the above parameter distance by $\|\mu_{i}-\hat{\mu}_{\pi(i)}\|+\|\Sigma_{i}-\hat{\Sigma}_{\pi(i)}\|_{F}\leq\epsilon,$ where $\|_{F}$ denotes the Frobenius norm.

Notationally, we will write $(\hat{F},\pi)\in\mathcal{D}_{\epsilon}(F)$ as shorthand for the statement that $\hat{F}$ is an $\epsilon$-correct subdivision for $F$ and $\pi$ is the (onto) function from $k$ to $k^{\prime}$ that groups $F$ into $\hat{F}$ as above.

Note that this definition, unlike the definition for $\epsilon$-close estimate, uses parameter distance as opposed to statistical distance. This is critical because our univariate algorithm will only be able to return an estimate that is an $\epsilon$-correct subdivision when the notion of “close” is in parameter distance, and not statistical distance because in general there could be a component of the univariate mixture of arbitrarily small variance, and we will only be able to match this to an additive guarantee and this implies nothing about the statistical distance between our estimate and the actual component.

A Robust Univariate Algorithm

In this section, we give a learning algorithm for univariate mixtures of Gaussians that will be the building block for our learning algorithm in $n$-dimensions. Unlike in the case of , our univariate algorithm will not necessarily be given a mixture of Gaussians for which all pairwise parameter distances are reasonably large. Instead, it could happen that we are given a mixture of (say) three Gaussians so that some pair has arbitrarily small parameter distance.

In the case in which we are guaranteed that all pairwise parameter distances are reasonably large, we can iterate the technical ideas in to give an inductive proof that a simple brute force search algorithm will return good estimates. We call this algorithm the Basic Univariate Algorithm. From this, we build a General Univariate Algorithm that will return a good estimate regardless of the parameter distances, although in order to do so we will need to relax the notion of a good estimate to something weaker: the algorithm return an $\epsilon$-correct subdivision.

In this section, we show that we can efficiently learn the parameters of univariate mixtures of Gaussians, provided that the components of the mixture have nonnegligible pairwise parameter distances. We refer to this algorithm as the Basic Univariate Algorithm. Such an algorithm will follow easily from Theorem 4—the polynomially robust identifiability of univariate mixtures. Throughout this section we will consider two univariate mixtures of Gaussians:

We will call the pair $F,F^{\prime}$ $\epsilon$-standard if $\sigma_{i}^{2},\sigma_{i}^{\prime 2}\leq 1$ and if $\epsilon$ satisfies:

$|\mu_{i}|,|\mu^{\prime}_{i}|\leq\frac{1}{\epsilon}$

$|\mu_{i}-\mu_{j}|+|\sigma_{i}^{2}-\sigma_{j}^{2}|\geq\epsilon$ and $|\mu^{\prime}_{i}-\mu^{\prime}_{j}|+|\sigma_{i}^{\prime 2}-\sigma_{j}^{\prime 2}|\geq\epsilon$ for all $i\neq j$

$\epsilon\leq\min_{\pi}\sum_{i}\left(|w_{i}-w_{\pi(i)}^{\prime}|+|\mu_{i}-\mu_{\pi(i)}^{\prime}|+|\sigma_{i}^{2}-\sigma_{\pi(i)}^{\prime 2}|\right)$, where the minimization is taken over all mappings $\pi:\{1,\ldots,n\}\rightarrow\{1,\ldots,k\}.$

There is a constant $c>0$ such that, for any $\epsilon$-standard $F,F^{\prime}$ and any $\epsilon<c$,

While the dependency on $k$ in Theorem 4 is very bad, as we show in Section 6, this exponential dependency on $k$ is necessary. Specifically, we give a construction of two $1/k$-standard distributions whose statistical distance is $O(e^{-k}).$

Given the polynomially robust identifiability guaranteed by the above theorem, and simple concentration bounds on the $i^{th}$ sample moment, it is easy to see that a brute-force search over a set of candidate parameter sets will yield an efficient algorithm that recovers the parameters for a univariate mixtures of Gaussians whose components have pairwise parameter distance at least $\epsilon$: roughly, the Basic Univariate Algorithm will take a polynomial number of samples, compute the first $4k-2$ sample moments, and compare those with the first $4k-2$ moments of each of the candidate parameter sets. The algorithm then returns the parameter set whose moments most closely match the sample moments. Theorem 4 guarantees that if the first $4k-2$ sample moments closely match those of the chosen parameter set, then the parameter set must be nearly accurate. To conclude the proof, we argue that a polynomial-sized set of candidate parameters suffices to guarantee that at least one set of parameters will yield moments sufficiently close to the sample moments, which, by simple concentration bounds, will be close to the true moments of the GMM. We state the corollary below, and defer the details of the algorithm and the proof of its correctness to Appendix C.

Suppose we are given access to independent samples from a GMM $\sum_{i=1}^{k}w_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2},x)$ with mean 0 and variance in the interval $[1/2,2],$ where $w_{i}\geq\epsilon$, and $|\mu_{i}-\mu_{j}|+|\sigma_{i}^{2}-\sigma_{j}^{2}|\geq\epsilon$. There exists a Basic Univariate Algorithm that, for any fixed $k$, has runtime at most $poly(\frac{1}{\epsilon},\frac{1}{\delta})$ samples and with probability at least $1-\delta$ will output mixture parameters $\hat{w}_{i},\hat{\mu}_{i},\hat{\sigma_{i}}^{2}$, so that there is a permutation $\pi:[k]\rightarrow[k]$ and

2 The General Univariate Algorithm

In this section we seek to extend the Basic Univariate Algorithm of Corollary 5 to the general setting of a univariate mixture of $k$ Gaussians without any requirements that the components have significant pair-wise parameter distances. In particular, given some target accuracy $\epsilon,$ and access to independent samples from a mixture $F$ of $k$ univariate Gaussians, we want to efficiently compute a mixture $F^{\prime}$ of $k^{\prime}\leq k$ Gaussians that is an $\epsilon$-correct subdivision of $F.$

There is a General Univariate Algorithm which, given $\epsilon,\delta>0$, and access to a GMM of $k$ Gaussians, $F=\sum_{i}w_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2})$ that is in near isotropic position and satisfies $w_{i}\geq\epsilon$, will run in time polynomial in $1/\epsilon$ and $1/\delta,$ and will return with probability at least $1-\delta$ a GMM of $k^{\prime}\leq k$ Gaussians $\hat{F}$ that is an $\epsilon$-correct subdivision of $F$.

The critical insight in building up such a General Univariate Algorithm is that if two components are actually close enough (in statistical distance), then the Basic Univariate Algorithm could never tell these two components apart from a single (appropriately) chosen Gaussian, because this algorithm only requires a polynomial number of samples. So given a target precision $\epsilon_{1}$ for the Basic Univariate Algorithm, there is some window that describes whether or not the algorithm will work correctly. If all pairwise parameter distances are either sufficiently large or sufficiently small, then the Basic Univariate Algorithm will function as if it were given sample access to a mixture that actually does meet the requirements of the algorithm. So as long as no parameter distance falls inside a particular window (which characterizes whether or not the algorithm will behave properly), the algorithm will return a correct computation.

However, when there is some parameter distance that falls inside the Basic Univariate Algorithm’s window, we are not guaranteed that the Basic Univariate Algorithm will fail safely. The idea, then, is to use many disjoint windows (each of which corresponds to running the Basic Univariate Algorithm with some target precision). If we choose enough such windows, each pairwise parameter distance can only corrupt a single run of the Basic Univariate Algorithm so a majority of the computations will be correct. We will never know which computations resulted from cases when no parameter distance fell inside the corresponding window, but we will be able to define a notion of consensus among these different runs of the Basic Univariate Algorithm so that a majority of the runs will agree, and any run which agrees with some computation that was correct will also be close to correct.

We defer the algorithm and proof of correctness to Appendix D

Partition Pursuit

In this section we demonstrate how to use the General Univariate Algorithm to obtain good additive approximations in $n$-dimensions. Roughly, we will project the $n$-dimensional mixture $F$ onto many close-by directions, and run the General Univariate Algorithm on each projection. This is also how the algorithm in is able to recover good additive estimates in $n$-dimensions. However we will have to cope with the additional complication that our univariate algorithm (the General Univariate Algorithm) does not necessarily return an estimate that is a mixture of $k$ Gaussians.

We explain in detail how the algorithm in is able to obtain additive approximation guarantees in $n$-dimensions, building on a univariate algorithm for learning mixtures of two Gaussians. Let $\epsilon_{3}>>\epsilon_{2}>>\epsilon_{1}$. Given any $\epsilon$-statistically learnable mixture of two Gaussians in $n$-dimensions, with high probability, for a direction $r$ chosen uniformly at random the parameter distance between the two Gaussians in $P_{r}[F]$ will be at least $\epsilon_{3}$. Then given such a direction $r$, we can choose $n^{2}$ different directions $r_{x,y}$ each of which are $\epsilon_{2}$-close to $r$ (i.e. $\|r-r_{x,y}\|\approx\epsilon_{2}$). The mean and variance of a component in $P_{u}[F]$ change continuously as we vary the direction $u$ from $r$ to $r_{x,y}$, and this implies that for $\epsilon_{2}<<\epsilon_{3}$, we will be able to consistently pair up estimates recovered from each projection, so that for each Gaussian we have $n^{2}$ different estimates in different directions of the projected mean and variance. Each of these estimates are accurate to within $\epsilon_{1}$ (i.e. this is the target precision that is given to the univariate algorithm). For any Gaussian, an estimate for the projected mean and the projected variance for a direction $r$ gives a linear constraint on the mean vector $\mu$ and the co-variance matrix $\Sigma$. As a result, if $\epsilon_{1}<<\epsilon_{2}$ then the precision is much finer than the condition number of this system of linear constraints on $\mu,\Sigma$ and this yields an accurate estimate in $n$-dimensions.

Let $\epsilon_{2},\epsilon_{1}>0$. Suppose $|m^{0}-\mu\cdot r|$,$|{m}^{ij}-\mu\cdot r^{ij}|$, $|v^{0}-r^{T}\Sigma r|$,$|v^{ij}-(r^{ij})^{T}\Sigma r^{ij}|$ are all at most $\epsilon_{1}$. Then Solve outputs $\hat{\mu}\in{\bf R}^{n}$ and $\hat{\Sigma}\in{\bf R}^{n\times n}$ such that $\|\hat{\mu}-\mu\|<\frac{\epsilon_{1}\sqrt{n}}{\epsilon_{2}}$, and $\|\hat{\Sigma}-\Sigma\|_{F}\leq\frac{6n\epsilon_{1}}{\epsilon_{2}^{2}}$. Furthermore, $\hat{\Sigma}\succeq 0$ and $\hat{\Sigma}$ is symmetric.

The algorithm to which this lemma refers is given in Appendix F.2

However, the General Univariate Algorithm does not always return a mixture of $k$ Gaussians, and can in fact return a mixture $\hat{F}^{u}$ of $k^{\prime}<k$ Gaussians provided that this mixture is still an $\epsilon_{1}$-correct subdivision of $P_{u}[F]$ (for some direction $u$). But then what happens if we consider two close-by directions, $u$ and $v$ and the number of Gaussians in the estimate $\hat{F}^{u}$ is different from the number of Gaussians in the estimate $\hat{F}^{v}$?

The key insight is that if we choose some direction $r$, and close-by directions $r_{x,y}$, if any estimate returned for $r_{x,y}$ has more components than the estimate returned for the direction $r$, then we have made progress because we have identified another Gaussian in the original mixture $F$. So here, rather than trying to use this estimate for $r_{x,y}$, we just start the algorithm over using $r_{x,y}$ as the original direction, and considering $n^{2}$ close-by directions.

The additional complication is that we must make sure every time we see a different number of components, that we’ve made progress. We can do so by maintaining a Window from $\epsilon_{1}$ to $\epsilon_{3}$, and we say that a Window is satisfied if the estimate $\hat{F}^{r}$ returned for some direction $r$ has all pairs of Gaussians either at parameter distance at least $\epsilon_{3}$, or at most the precision $\epsilon_{1}$ of the General Univariate Algorithm. Then if we consider close-by directions $r_{x,y}$ (that are $\epsilon_{2}$-close to $r$, for $\epsilon_{1}<<\epsilon_{2}<<\epsilon_{3}$), we can ensure that whenever we see a different number of components in the estimate corresponding to some direction $r_{x,y}$, there are more components. When we see more components, we may need to shift the Window $W$ to a Window $W^{\prime}$ so that in this new direction $r_{x,y}$, the Window $W^{\prime}$ is satisfied. We take $r_{x,y}$ as the new base direction. But we have made progress because we have identified a new component in the mixture.

We state our main theorem in this section, and defer the algorithm and proof to Appendix F

Given an $\epsilon$-statistically learnable GMM $F$ in isotropic position, the Partition Pursuit Algorithm will recover an $\epsilon$-correct sub-division $\hat{F}$ and if $F$ has more than one component, $\hat{F}$ also has more than one component.

Clustering and Recursion

In this section, we give an efficient algorithm for learning an estimate $\hat{F}$ that is $\epsilon$-close to the actual mixture $F$. Partition Pursuit assumes that the mixture $F$ is in isotropic position, and even though $F$ is not necessarily in isotropic position, we will be able to get around this hurdle by first taking enough samples to compute a transformation that places the mixture $F$ in nearly isotropic position and then applying this transformation to each sample from the oracle. The main technical challenge in this section is actually what to do when the mixture $\hat{F}$ returned by Partition Pursuit is a good additive approximation to $F$ (i.e. it is an $\epsilon_{1}$-correct subdivision with $\epsilon_{1}<<\epsilon$), but is not $\epsilon$-close to the mixture $F$. This can only happen if there is a component in $F$ that has a very small variance in some direction. Consider for example, two Gaussians in one dimension $\mathcal{N}(0,\gamma)$ and $\mathcal{N}(0,\gamma+\epsilon_{1})$. Even if $\epsilon_{1}$ is very small, if $\gamma$ is much smaller, then the statistical distance between these two Gaussians can be arbitrarily close to $1$.

So the high-level idea is that if the estimate $\hat{F}$ returned by Partition Pursuit is not $\epsilon$-close to $F$ (but $\hat{F}$ is an $\epsilon_{1}$-correct subdivision of $F$ for $\epsilon_{1}<<\epsilon$), then it must be the case that some component $F_{i}$ of $F$ has a co-variance matrix $\Sigma_{i}$ so that for some direction $v$, $v^{T}\Sigma_{i}v$ is very small. Then we can use this direction $v$ to still make progress: If we project the mixture $F$ onto $v$, we will be able to cluster accurately. There will be some partition of the Gaussians in $F$ into two disjoint, non-empty sets of components $S,T$ and some clustering scheme that can accurately clusters points sampled from $F$ into points that originated from a component in $S$ and points that originated from a component in $T$. So we can hope to accurately cluster enough points sampled from $F$ into sets of points that originated from $S$ and sets of points that originated from $T$, and then we can run our learning algorithm (with a smaller maximum of at most $k-1$ components) on each set of points. By induction, this learning algorithm will return close estimates, and if we take a convex combination of these estimates we obtain a new estimate $\hat{F}^{\prime}$ that is $\epsilon$-close to $F$. The main technical challenge is in showing that if there is some component of $F$ with a small enough variance in some direction $v$, then we can accurately cluster points sampled from $F$. Given this, our main result follows almost immediately from an inductive argument.

2 How to Cluster

Here we give formalize the notion of a clustering scheme. Additionally, we state the key lemmas that will be useful in showing that if $\hat{F}$ is not an $\epsilon$-close estimate to $F$, then we can use $\hat{F}$ to construct a good clustering scheme that makes progress on our learning problem.

We will call $A,B\subset\Re^{n}$ a clustering scheme if $A\cap B=\emptyset$

For $A\subset\Re^{n}$, we will write $P[F_{i},A]$ to denote $Pr_{x\sim F_{i}}[x\in A]$ - i.e. the probability that a randomly chosen sample from $F_{i}$ is in the set $A$.

If we have a direction $v$ and some component $\hat{F}_{i}$ which has small variance in direction $v$, we want to use this direction to cluster accurately. The intuition is clearest in the case of mixtures of two Gaussians: Suppose one of the components, say $\hat{F}_{1}$, had small variance on direction $v$. If the entire mixture is in isotropic position, then the variance of the mixture when projected onto direction $v$ is $1$. This can only happen if either the difference in projected means $|v^{T}(\hat{\mu}_{1}-\hat{\mu}_{2})|$ is large or the variance of $\hat{F}_{2}$ on direction $v$ is large. In the first case, we can choose an interval around each projected (estimate) mean $v^{T}\hat{\mu}_{1}$ and $v^{T}\hat{\mu}_{2}$ so that with high probability, any point sampled from $F_{1}$ is contained in the interval around $v^{T}\hat{\mu}_{1}$ and similarly for $F_{2}$. If, instead, the variance of $F_{2}$ when projected onto $v$ is large, then again a small interval around the point $v^{T}\hat{\mu}_{1}$ will contain most samples from $F_{1}$, but because the maximum density of $v^{T}F_{2}$ is never large and the interval around $v^{T}\hat{\mu}_{1}$ is not too large either, most samples from $F_{2}$ will not be contained in the interval. This idea is the basis of our clustering lemmas, although there will be additional complications when the mixture contains more than two Gaussians, the intuition is close to the same.

Let $(\hat{F},\pi)\in\mathcal{D}_{\epsilon_{1}}(F)$. Suppose also that $\hat{F}$ is a mixture of $k^{\prime}$ components.

Suppose that for some direction $v$, for all $i$: $v^{T}\hat{\Sigma}_{i}v\leq\epsilon_{2}$, for $\epsilon_{1}\leq\frac{\sqrt{\epsilon_{2}}}{2\epsilon_{3}}$. If there is some bi-partition $S\subset[k^{\prime}]$ s.t. $\forall_{i\in S,j\in[k^{\prime}]-S}|v^{T}\hat{\mu}_{i}-v^{T}\hat{\mu}_{j}|\geq\frac{3\sqrt{\epsilon_{2}}}{\epsilon_{3}}$ then there is a clustering scheme $(A,B)$ (based only on $\hat{F}$) so that for all $i\in S,j\in\pi^{-1}(i)$, $P[F_{i},A]\geq 1-\epsilon_{3}$ and for all $i\notin S,j\in\pi^{-1}(i)$, $Pr[F_{i},B]\geq 1-\epsilon_{3}$.

This lemma corresponds to the first case in the above thought exercise when there is some bi-partition of the components so that all pairs of projected means across the bi-partition are reasonably separated.

Suppose that for some direction $v$ and some $i\in[k^{\prime}]$ such that: $v^{T}\hat{\Sigma}_{i}v\leq\epsilon_{m}$, for $\epsilon_{m}>>\epsilon_{1}$. If there is some bi-partition $S\subset[k^{\prime}]$ s.t.

(and $\epsilon_{t}<<\epsilon_{3}^{3}$) then there is a clustering scheme $A,B$ such that for all $i\in S,j\in\pi^{-1}(i)$, $P[F_{i},A]\geq 1-\epsilon_{3}$ and for all $i\notin S,j\in\pi^{-1}(i)$, $Pr[F_{i},B]\geq 1-\epsilon_{3}$.

This lemma corresponds to the second case to the second case, when there is some bi-partition of the components so that one side of the bi-partition has projected variances that are much larger than the other.

The proofs of these lemmas, along with additional technical details are given in Appendix G.2

3 Making Progress when there is a Small Variance

We state a lemma from which formalizes the intuition that if there is no component in $\hat{F}$ with small variance in any direction, the $\hat{F}$ is a good statistical estimate to $F$:

Suppose $\|\hat{\mu}_{i}-\mu_{i}\|\leq\epsilon_{1}$, $\|\hat{\Sigma}_{i}-\Sigma_{i}\|_{F}\leq\epsilon_{1}$, and $|\hat{w}_{i}-w_{i}|\leq\epsilon_{1}$, if either $\|\Sigma^{-1}_{i}\|_{2}\leq\frac{1}{2\epsilon_{m}}$ or $\|\hat{\Sigma}^{-1}_{i}\|_{2}\leq\frac{1}{2\epsilon_{m}}$ then

We will use this lemma as a building block to prove:

The Hierarchical Clustering Algorithm either returns an $\epsilon$-close statistical estimate $\hat{F}$ for $F$, or returns a clustering scheme $A,B$ such that there is some bipartition $S\subset[k]$ such that for all $i\in S,j\in\pi^{-1}(i)$, $P[F_{i},A]\geq 1-\epsilon_{3}$ and for all $i\notin S,j\in\pi^{-1}(i)$, $Pr[F_{i},B]\geq 1-\epsilon_{3}$. And also $S,[k]-S$ are both non-emtpy.

We defer the algorithm and the proof of correctness to Appendix G.

4 Recursion

[Isotropic Projection Lemma] Given a mixture of $k$ $n$-Dimensional Gaussians $F=\sum_{i}w_{i}F_{i}$ that is in isotropic position and is $\epsilon$-statistically learnable, with probability $\geq 1-\delta$ over a randomly chosen direction $u$, there is some pair of Gaussians $F_{i},F_{j}$ s.t. $D_{p}(P_{u}[F_{i}],P_{u}[F_{j}])\geq\frac{\epsilon^{5}\delta^{2}}{50n^{2}}$.

We defer a proof of this lemma to Appendix H

Let $H_{a}(\epsilon,\delta,k),H_{i}(\epsilon,\delta,k)$ be the inverse of the number of samples needed by the High Dimensional Anisotropic Algorithm and the High Dimensional Isotropic Algorithm respectively when given target precision $\epsilon$ (and access to an $\epsilon$-statisically learnable distribution), an upper bound $k$ on the number of Gaussians, and an error parameter $\delta$.

Given $k,\epsilon$, and a mixture of at most $k$ Gaussians $F$ that is $\epsilon$-statistically learnable High Dimensional Anisotropic Algorithm returns an estimate $\hat{F}$ that is $\epsilon$-close to the actual mixture $F$.

We prove this theorem by induction. Let $\epsilon_{k-1}=H_{a}(\frac{\epsilon}{2},\delta,k-1)$.

We assume by induction that both the High Dimensional Isotropic Algorithm and the High Dimensional Anisotropic Algorithm return an $\epsilon$-close estimate for all values of $k^{\prime}\leq k-1$. We then consider both algorithms for the case of $k$:

Consider the High Dimensional Isotropic Algorithm which is given $k,\epsilon$, and a mixture of at most $k$ Gaussians $F$ that is $\epsilon$-statistically learnable and is in isotropic position: We first run the Hierarchical Clustering Algorithm with parameters $\epsilon,\delta,\epsilon_{3},k$ where $\epsilon_{3}=\frac{1}{2}\epsilon\epsilon_{k-1}\delta$. If this algorithm returns an estimate $\hat{F}$, we can return this estimate and it is guaranteed to be $\epsilon$-close to the actual mixture.

Note that if the number of components in $F$ is $1$, then the Hierarchical Clustering Algorithm will necessarily return an estimate $\hat{F}$, because there is no partitioning scheme that partitions $F$ into two subsets of components that are both non-empty. This establishes the base case in the inductive argument.

Otherwise the output of the Hierarchical Clustering Algorithm is a clustering scheme $(A,B)$ with the property that there is some partition $S,T$ of the Gaussians in $F$ ($S,T\neq\emptyset$) and for all $i\in S$, $Pr_{x\sim F_{i}}[x_{i}\in A]\geq 1-\epsilon_{3}$, and $j\in T$, $Pr_{x\sim F_{i}}[x_{i}\in B]\geq 1-\epsilon_{3}$. Let $F_{S},F_{T}$ be the (re-weighted) mixtures that result from placing every component in $S$ from $F$ into $F_{S}$, and every component in $T$ from $F$ into $F_{T}$. Note that $F_{S},F_{T}$ are still $\epsilon$-statistically learnable, but may not be in isotropic position any longer.

All samples $x_{1},x_{2},...,x_{m}$ are either in $A$ or $B$

The number of samples in $A$ and the number of samples in $B$ will each be at least $\frac{1}{\epsilon_{k-1}}$

All samples are clustered correctly - i.e. if $x_{i}\in A$, then $x_{i}$ was generated by some Gaussian $F_{j}$ with $j\in S$ and if $x_{i}\in B$, then $x_{i}$ was generated some Gaussian $F_{j}$ with $j\in T$.

Let $X_{S},X_{T}$ be the samples from $x_{1},x_{2},...,x_{m}$ that are in $A,B$ respectively. We can then run the the High Dimensional Anisotropic Algorithm with parameters $\frac{\epsilon}{2},\delta,k-1$ on each set $X_{S}$ and $X_{T}$. Let the algorithm return the mixtures $\hat{F}_{A},\hat{F}_{B}$ respectively. We return a convex combination of these mixtures, $\hat{F}=\frac{|X_{S}|}{m}\hat{F}_{A}+\frac{|X_{T}|}{m}\hat{F}_{B}$. The estimates $\hat{F}_{A},\hat{F}_{B}$ are $\epsilon$-close estimates to $F_{S},F_{T}$ respectively. We can write $F=w_{A}F_{A}+w_{B}F_{B}$, and with high probability $\frac{|X_{S}|}{m}$, $\frac{|X_{T}|}{m}$ will be close to $w_{A},w_{B}$ respectively. Then this implies that $\hat{F}$ is $\epsilon$-close to $F$. Thus by induction, the output of the High Dimensional Isotropic Algorithm is an estimate $\hat{F}$ that is $\epsilon$-close to $F$.

We need to also verify by induction that the output of the High Dimensional Anisotropic Algorithm is also an $\epsilon$-close estimate to $F$. So suppose that the input to the High Dimensional Anisotropic Algorithm is a mixture of at most $k$ Gaussians, that is $\epsilon$-statistically learnable and is not necessarily in isotropic position.

Exponential Dependence on k𝑘k is Inevitable

In this section, we present a lower bound, showing that the inverse exponential dependency on the number of Gaussian components in each mixture is necessary, even for mixtures in just one dimension. We show this by giving a simple construction of two 1-dimensional distributions, $D_{1},D_{2}$ that are $1/(2m)$-standard. Specifically, they are mixtures of at most $m$ Gaussians, such that the weights of all components of each mixture are at least $1/(2m)$, and the parameter distance between the pair of distributions is at least $1/(2m),$ but $||D_{1}-D_{2}||_{1}\leq e^{-m/30},$ for sufficiently large $m$. The construction hinges on the inverse exponential (in $k\approx\sqrt{m}$) statistical distance between $\mathcal{N}(0,2),$ and the mixtures of infinitely many Gaussians of unit variance whose components are centered at multiples of $1/k$, with the weight assigned to the component centered at $i/k$ being given by $N(0,1,i/k).$ Verifying that this is true is a straight-forward exercise in Fourier analysis. The final construction truncates the mixture of infinitely many Gaussians by removing all the components with centers a distance greater than $k$ from 0. This truncation clearly has negligibly small effect on the distribution. Finally, we alter the pair of distributions by adding to both distributions, Gaussian components of equal weight with centers at $-k^{2}/k,(-k^{2}+1)/k,\ldots,k,$ which ensures that in the final pair of distributions, all components have significant weight.

There exists a pair $D_{1},D_{2}$ of $1/(4k^{2}+2)$-standard distributions that are each mixtures of $k^{2}+1$ Gaussians such that

We can define $F_{k}(x)^{N}=c_{k}\sum_{i=-N}^{N}\frac{1}{\sqrt{\pi}}e^{-(i/k)^{2}}\mathcal{N}(i/k,1/2,x),$, and we give a plot of $F_{k}^{N}$ for $k=2,N=2$ in Figure 1a and the corresponding plot of each component, and in Figure 1b we give a plot of each component of $F_{k}^{N}$ for $k=4,N=8$.

Conclusions

We give an estimator that converges to the true distribution at an inverse polynomial rate, and this result has implications for polynomial-time clustering and density estimation. A natural question is: “What is the optimal rate of convergence?” This question is wide open, and all we can say for certain is that the rate of convergence is at worst polynomial in the dimension and the inverse of the desired accuracy, and exponential in the number of components. We made no attempt here to optimize the constants in the exponent of the rate of convergence and even if we had, the theoretical runtime would still be extremely impractical. This, however, raises the practically relevant question of whether aspects of our algorithm can be combined with existing heuristics that seem to perform well in most applications. For example, the brute-force-search component of our univariate algorithm is clearly expensive; perhaps employing existing heuristics (such as the EM algorithm) for the univariate problems, in conjunction with aspects of our dimension-reduction machinery might yield improved efficiency on real-world instances.

Additionally, we note that much of the machinery we developed—from the “deconvolution” argument for the polynomially robust identifiability, to the partition pursuit and hierarchical clustering for the dimension reduction arguments, seem to be relatively general and robust. We suspect that such tools could be applied to yield corresponding results for other families of distributions.

Acknowledgements

We are grateful to Paul Valiant, for suggesting the lower-bound construction of Section 6 and many helpful discussions throughout; and are indebted to Adam Tauman Kalai for introducing us to this beautiful line of research, and for all his guidance, encouragement and deep insights about mixtures of Gaussians.

References

Appendix A In-Depth Outline

To start, suppose that we are given access to independent samples from a mixture of $k$ Gaussians, and given a unit vector $r$ with the following promise: for each pair $G_{1},G_{2}$ of components, in the projection of the mixture onto $r$, either the projections of $G_{1}$ and $G_{2}$ have reasonably different parameters ($>\epsilon_{2}$), or their projections are so close that our algorithm could never tell them apart from a single Gaussian (parameter distance at most $\epsilon_{0}<<\epsilon_{1},$ where $\epsilon_{1}<<\epsilon_{2}$ is the desired accuracy of the 1-d parameter learning algorithm. In this case, our 1-d parameter recovery algorithm will perform correctly, and return some $\epsilon_{1}$-accurate parameters for a mixture of $k^{\prime}\leq k$ components.

Thus in general, for a given desired accuracy $\epsilon_{1}$, there is some critical window, namely $[\epsilon_{0},\epsilon_{2}],$ associated with the 1-d learning algorithm that determines if it will function correctly. In a given projection, as long as no pair of components have parameter distances that fall within this window, then any pair of Gaussians is either reasonably different in parameters, or so close in parameters that the algorithm will never be able to tell the difference.

In this way, if an algorithm designer is told the parameters of a given mixture of $k$ Gaussians, he could construct an algorithm that would have been able to find some of these parameters. The algorithm would project onto a random direction $r$, and based on the pairwise parameter distances between the univariate Gaussians, there will be some window (i.e. some choice of a target precision with which to run the algorithm), bounded below by some polynomial in the desired output accuracy, so that the algorithm would function correctly. The problem is that while there is always some window that would work for any mixture of $k$ univariate Gaussians, we don’t know what window to use, and in general if we run the algorithm on a bad window, we aren’t guaranteed that the algorithm will fail in a safe way.

To get around this, we run the 1-d Learning Algorithm algorithm many times on different windows that do not intersect. Because there are only $k$ univariate Gaussians, and thus at most $k^{2}$ different distances between component parameters in any given projection, at most $k^{2}$ of these windows can be corrupted. If we choose sufficiently many windows (but still a polynomial number), a majority of the windows will yield correct parameters. Even though we can never determine which windows were good and which were bad, we can return the parameters generated by some window in consensus with a majority, and in this way, regardless of whether the window was good or bad, it is in consensus with a good window and must also be close to the correct parameters.

It is important to stress that even after the above consensus is conducted on a given projection, we still cannot be guaranteed that our univariate algorithm returns a mixture of $k$ Gaussians. Instead, it will return some mixture of $k^{\prime}<k$ Gaussians, where an element in the mixture might correspond to (say) a pair of Gaussians in the original mixture that were too close to differentiate in the given projection.

A.2 Partition Pursuit

This brings us to the second obstacle outlined in Section 1.4: in order to recover the $n$-dimensional parameters, we will need estimates of the parameters of the Gaussians when projected on many different directions. But, as mentioned above, the univariate algorithm will not necessarily return a mixture of $k$ Gaussians, and even if we choose a direction $r^{\prime}$ that is sufficiently close to $r$ (but still $|r-r^{\prime}|>>\epsilon_{1},$ the accuracy of the 1-d algorithm), it may be the case that the univariate algorithm for direction $r$ returned $k^{\prime}$ Gaussians, and the univariate algorithm for direction $r$ returned $k^{\prime\prime}\neq k^{\prime}$ Gaussians. How do we pair up these estimates in a consistent manner?

The key insight is that we are actually making progress if we see more Gaussians when projecting onto a different direction. If we choose a new direction, and we see a mixture of Gaussians with more components, we should backtrack and start over as if this was the direction we originally chose. We may have to slide the window corresponding to our 1-d algorithm and learn at a finer precision than what we chose previously, but this finer precision will still be polynomially bounded. Effectively, we are clustering the Gaussian components into clusters with the property that the components of each cluster are indistinguishable in each of the one-dimensional projections that we have considered. In order to make this idea work properly, we also need to ensure that we maintain a minimum parameter distance between all Gaussians clusters that we have seen (i.e. this distance is much larger than our 1-d accuracy $\epsilon_{1}$), so that when we choose a new direction $r^{\prime}$ sufficiently close to $r,$ Gaussian component cannot switch clusters. Thus at each stage, each cluster of Gaussians either continues to be a cluster, or it gets partitioned into several clusters of Gaussians.

A.3 Hierarchical Clustering

The final obstacle outlined in Section 1.4 can be addressed easily via an accurate clustering of the input samples together with a $k$-Gaussian analog of the Random Projection Lemma. Intuitively, the only way that a set of high-dimensional Gaussians with significant statistical distance, when projected onto a random vector, will appear nearly identical is if the re-weighted mixture of the Gaussians in this set is very far from isotropic position. This motivates the hope that if we have recovered some mixture of $k^{\prime}<k$ components, then it must be the case that whichever of these components contains multiple original Gaussians has covariance matrix very far from isotropic. Thus such a component must have at least one very small eigenvalue. Given the eigenvector corresponding to such an eigenvalue, we should be able to very accurately cluster the sample points into some partition of the original Gaussians. This motivates the following slightly more specific version of the high-level algorithm approach: Given that we have recovered parameters for a mixture of $k^{\prime}<k$ components:

Learn the parameters of some mixture of $k^{\prime}\leq k$ Gaussians, where each learned Gaussian component corresponds to one or more of the Gaussians in the original mixture.

If $k^{\prime}<k$, for each of the $k^{\prime}$ components recovered in the previous step:

If the $i^{th}$ component has covariance matrix “not too far” from isotropic, then conclude that it corresponds to a single Gaussian in the original mixture.

there is a very small eigenvalue of the covariance matrix, so project the sample points onto the corresponding eigenvector, and accurately cluster the sample points that come from this component

Given the sample points corresponding to one of the components, rescale these data points so this component (which was very far from isotropic), is now in isotropic position, and repeat the entire algorithm on this sub-mixture

The final observation that guarantees that our algorithm will make progress with every iteration, and thus terminate after a polynomial number of steps is the following analog of the Random Projection Lemma for the $k$-Gaussians setting. Given a mixture of $k$ Gaussians in isotropic position, with high probability over random unit vectors $r$, there will be some pair of projected Gaussians whose parameters are reasonably different. Thus, in every projection, we will, with high probability, see what appears to be a mixture of at least two components.

Appendix B Polynomially Robust Identifiability

We now sketch the rough outline of the proof of Theorem 4. While there are considerable technical details, the main proof ideas are identical to those used in to prove the analogous theorem in the case that $n=k=2.$

Our proof will be via induction on $\max(n,k)$. We start by considering the constituent Gaussian of minimal variance in the mixtures. Assume without loss of generality that this minimum variance component is the first component of $F,$ and denote it by $N_{1}$. If there is no component of $F^{\prime}$ whose mean, variance, and mixing weight very closely matches those of $N_{1}$, then we argue that there is a significant disparity in the low order moments of $F$ and $F^{\prime}$, no matter what the other Gaussian components are. (This argument is rather involved, and we will give the high-level sketch in the next paragraph.) If there is a component $N_{1}^{\prime}$ of $F^{\prime}$ whose mean, variance, and mixture weight very closely matches those of $N_{1}$, then we argue that we can remove $N_{1}$ from $F$ and $N_{1}^{\prime}$ from $F^{\prime}$ with only negligible effect on the discrepancy in the low-order moments. More formally, let $H$ be the mixture of $n-1$ Gaussians obtained by removing $N_{1}$ from $F$, and rescaling the weights so as to sum to one, and define $H^{\prime},$ a mixture of $k-1$ Gaussians analogously. Then, assuming that $N$ and $N^{\prime}$ are very similar, the disparity in the low-order moments of $H$ and $H^{\prime}$ is almost the same as the disparity in low-order moments of $F$ and$F^{\prime}$. We can then apply the induction hypothesis to the mixtures $H$ and $H^{\prime}$.

We now return to the problem of showing that if the skinniest Gaussian in $F$ cannot be paired with a component of $F^{\prime}$ with similar mean, variance, and weight, that there must be a polynomially-significant discrepancy in the low-order moments of $F$ and $F^{\prime}$. This step relies on ’deconvolving’ by a Gaussian with an appropriately chosen variance (this corresponds to running the heat equation in reverse for a suitable amount of time). We define the operation of deconvolving by a Gaussian of variance $\alpha$ as $\mathcal{F}_{\alpha}$; applying this operator to a mixture of Gaussians has a particularly simple effect: subtract $\alpha$ from the variance of each Gaussian in the mixture (assuming that each constituent Gaussian has variance at least $\alpha$). If $\alpha$ is negative, this is just convolution.

Let $F(x)=\sum_{i=1}^{n}w_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2},x)$ be the probability density function of a mixture of Gaussian distributions, and for any $\alpha<\min_{i}\sigma_{i}^{2},$ define

The key step will be to show that if the skinniest Gaussian in either of the two mixtures cannot be paired with a nearly identical Gaussian in the other mixture, then there is some $\alpha$ for which the resulting mixtures, after applying the operation $\mathcal{F}_{\alpha}$, have large statistical distance. Intuitively, this deconvolution operation allows us to isolate Gaussians in each mixture and then we can reason about the statistical distance between the two mixtures locally, without worrying about the other Gaussians in the mixture.

Given this statistical distance between the transformed pair of mixtures, we the fact that there are relatively few zero-crossings in the difference in probability density functions of two mixtures of Gaussians (Proposition 19) to show that this statistical distance gives rise to a discrepancy in at least one of the low-order moments of the pair of transformed distributions. To complete the argument, we then show that applying this transform to a pair of distributions, while certainly not preserving statistical distance, roughly preserves the combined disparity between the low-order moments of the pair of distributions. The complete proof can be found in Appendix B.

B.2 Theorem 4

In this section we give the complete proof of the polynomially robust identifiability of univariate mixtures of $k$ Gaussians (Theorem 4). For convenience, we restate the theorem and all necessary definitions. We make frequent reference to the simple properties of Gaussians and tail bounds provided in Appendix J. Throughout this section we will consider two univariate mixtures of Gaussians:

Definition 6. We will call the pair $F,F^{\prime}$ $\epsilon$-standard if $\sigma_{i}^{2},\sigma_{i}^{\prime 2}\leq 1$ and if $\epsilon$ satisfies:

$|\mu_{i}|,|\mu^{\prime}_{i}|\leq\frac{1}{\epsilon}$

$|\mu_{i}-\mu_{j}|+|\sigma_{i}^{2}-\sigma_{j}^{2}|\geq\epsilon$ and $|\mu^{\prime}_{i}-\mu^{\prime}_{j}|+|\sigma_{i}^{\prime 2}-\sigma_{j}^{\prime 2}|\geq\epsilon$ for all $i\neq j$

$\epsilon\leq\min_{\pi}\sum_{i}\left(|w_{i}-w_{\pi(i)}^{\prime}|+|\mu_{i}-\mu_{\pi(i)}^{\prime}|+|\sigma_{i}^{2}-\sigma_{\pi(i)}^{\prime 2}|\right)$, where the minimization is taken over all mappings $\pi:\{1,\ldots,n\}\rightarrow\{1,\ldots,k\}.$

Theorem 4. There is a constant $c>0$ such that, for any $\epsilon$-standard $F,F^{\prime}$ and any $\epsilon<c$,

The following definition of the deconvolution operation will be central to our proof of Theorem 4:

Definition 10. Let $F(x)=\sum_{i=1}^{n}w_{i}\mathcal{N}(\mu_{i},\sigma_{i}^{2},x)$ be the probability density function of a mixture of Gaussian distributions, and for any $\alpha<\min_{i}\sigma_{i}^{2},$ define

The following lemma argues that if the skinniest Gaussian in mixture $F$ can not be matched with a sufficiently similar component in the mixture $F^{\prime}$, then there is some $\alpha$, possibly negative, such that $\max_{x}|\mathcal{F}_{\alpha}(F)(x)-\mathcal{F}_{\alpha}(F^{\prime})(x)|$ is significant. Furthermore, every component in the transformed mixtures have variances that are not too small.

Suppose $F,F^{\prime}$ are $\epsilon$-standard. Suppose without loss of generality that the Gaussian of minimal variance is $\mathcal{N}(\mu_{1},\sigma_{1}^{2}),$ and there is some $\gamma$ satisfying $\epsilon/4>\gamma>0$ such that for all $i,$ at least one of the following holds:

$|\sigma_{1}^{2}-\sigma_{i}^{\prime 2}|>\gamma^{8}$

Then there is some $\alpha$ such that either

$\max_{x}(|\mathcal{F}_{\alpha}(F)(x)-\mathcal{F}_{\alpha}(F^{\prime})(x)|)\geq\frac{1}{2\gamma\sqrt{2\pi}}$ and the minimum variance in any component of $\mathcal{F}_{\alpha}(F)$ or $\mathcal{F}_{\alpha}(F^{\prime})$ is at least $\gamma^{4},$

$\max_{x}(|\mathcal{F}_{\alpha}(F)(x)-\mathcal{F}_{\alpha}(F^{\prime})(x)|)\geq\frac{1}{2\gamma^{8}\sqrt{2\pi}}$ and the minimum variance in any component of $\mathcal{F}_{\alpha}(F)$ or $\mathcal{F}_{\alpha}(F^{\prime})$ is at least $\gamma^{18}.$