Learning Mixtures of Spherical Gaussians via Fourier Analysis

Designing algorithms that estimate the parameters of an underlying probability distributions is a central theme in statistical learning theory. An important special instance of this learning problem is the case when the underlying distribution is known to be a finite mixture of Gaussian distributions in $d$-dimensional euclidean space, because such mixtures are popular models for high-dimensional data clustering and learning mixture of Gaussians in an unsupervised setting has been a topic of intensive research for the last few decades.

In its most explicit form, the underlying problem is as follows: we have access to random samples drawn independently as per some Gaussian mixture $\mu:=\omega_{1}\mu_{1}+\cdots+\omega_{k_{0}}\mu_{k_{0}}$, where $(\omega_{1},\cdots,\omega_{k_{0}})$ is some probability vector with strictly positive components, and each $\mu_{l}$ is a Gaussian density in $\mathbb{R}^{d}$, with mean $y_{l}\in\mathbb{R}^{d}$ and covariance $\Sigma_{l}$. The main task is to estimate the parameter set $\{(\omega_{1},y_{1},\Sigma_{1}),\cdots,(\omega_{k_{0}},y_{k_{0}},\Sigma_{k_{0}})\}$ of the density function $\mu$, within a pre-specified accuracy $\delta$. For the purpose of this paper, we will restrict to the scenario where all the Gaussian components are spherical.

Many of the earlier approaches to this learning problem were based on local search heuristics, $e.g.$ the EM algorithm and $k$-means heuristics, that resulted in weak performance guarantees. In [6], Dasgupta presented the first provably correct algorithm for learning a mixture of Gaussians, with a common unknown covariance matrix, using only polynomially many (in dimension as well as number of components) samples and computation time, under the assumption that the minimum separation between the centers of component Gaussians is at least $\Omega(\mathrm{polylog}(kd)\sqrt{d})$. In a consequent work, Dasgupta and Schulman ([4]) showed a variation of EM to work with minimum separation only $\mathrm{polylog}(kd)d^{\frac{1}{4}}$, when the components are all spherical. Subsequently, many more algorithms (see [12] and the references therein) with improved sample complexity and/or computational complexity have since been devised, that work on somewhat weaker separation assumption; in particular, the SVD-based algorithm of Vempala and Wang [14] learns a mixture of spherical Gaussians with poly-sized samples and polynomially many steps, under the separation assumption

We note that when $k=\Theta(2^{d})$, the above separation requirement translates to minimum separation of $\Omega(d^{\frac{3}{4}})$. In another line of work, Kalai-Moitra-Valiant [10] and Moitra-Valiant [12], the question of polynomial learnability of arbitrarily separated mixture of Gaussians (and more general families of distributions have been investigated by Belkin and Sinha in [2]). It has been established that, for Gaussian mixture with a fixed number of components and the centers having a minimum required separation, there is an algorithm that runs in polynomial time and uses polynomially many samples to learn the parameters of the mixture to any desired accuracy, with arbitrarily high probability.

Recently, Hopkins and Li ([7]), and Diakonikolas et al ([5]) devised polynomial time learning algorithms that work with minimum separation of $k^{\epsilon}$ (for any $\epsilon>0$). In [13], Regev and Vijayaraghavan considered the question of obtaining a lower bound on the separation of centers necessary for the polynomial learnability of Gaussian mixtures. They devised an iterative algorithm for amplifying accuracy of parameter estimation that, given initializer $y_{1}^{0},\cdots,y_{k}^{0}$ and a desired accuracy parameter $\delta>0$, uses polynomially many samples from the underlying mixture and computation steps to return $y_{1}^{\prime},\cdots,y_{k}^{\prime}$, that lie within Hausdorff distance at most $\delta$ from the true centers; for more details, see Theorem 4.1 of the full version of [13]. One of their results establishes that, in constant dimension $d=O(1)$, with minimum separation at least $\Omega(1)$, any uniform mixture of spherical Gaussians can be learned to desired accuracy $\delta$ (with high probability) in number of samples and computation time depending polynomially on the number of components and $\delta^{-1}$.

In the present paper we consider the question of learning the centers of mixture of Gaussians, with minimum separation $\Omega(\sqrt{d})$, in number of samples and computational time that depends polynomially on the ambient dimension $d$, and the number of components $k$. Here is the main theorem:

Main Theorem: There is an absolute constant $c_{0}\in(0,1)$ such that, given a mixture of at most $k$ spherical Gaussians in $\mathbb{R}^{d}$, with known weights and known identical variance $\sigma^{2}$, for which all the mixing weights are in $[c,c^{-1}]$ for some arbitrary absolute constant $c\in(0,1)$, and the minimum separation of the centers of the components satisfies

there is a randomized algorithm that recovers (with high probability) the centers of the mixture up to accuracy $\delta$ in time $\mbox{poly}(k,d,\delta^{-1},\sigma^{-1})$.

In particular, we generalize the main upper bound in [13] (see theorem 5.1 of [13]) to arbitrary dimension ($d$) and number of components ($k$), when the variance and weights of the underlying mixture is known. We notice that our minimum separation requirement is independent of $k$, the number of components, and for $d=O(\log k)$, this is weak minimum separation requirement.

In what follows, we will assume that $\sigma=1$; by scale-invariance of the condition in the theorem above, this is without loss of generality.

We first observe that if two clusters of centers are very far (i. e. the minimum distance of a center in one cluster is at least $\sqrt{d}$ away from every center of the other cluster), then the samples are unambiguously of one cluster only. This is shown in Lemma 1, allowing us to reduce the question to the case when a certain proximity graph — defined on the centers — is connected.

In algorithm $LearnMixture$, we use the Johnson-Lindenstrauss lemma and project the data in the ambient space to $d+1$ carefully chosen subspaces of dimension at most $O(1+\min(\log k,d))$, and show that if we can separate the Gaussians in these subspaces, the resulting centers can be used to obtain a good estimate of the centers in the original question. Thus the question is further reduced to one where the dimension $d=O(\log k)$.

Next, in Lemma 2, we show that if the number of samples is chosen appropriately, then all the $k_{0}$ centers are with high probability contained in a union of balls $B_{i}$ of radius $2\sqrt{d}$ centered around the data points. This allows us to focus our attention to $n$ balls of radius $2\sqrt{d}$.

The main idea is that in low dimensions, it is possible to efficiently implement a deconvolution on a mixture of identical spherical Gaussians having standard deviation $\sigma=1$, and recover a good approximation to a mixture of gaussians with the same centers and weights, but in which each component now has standard deviation $\bar{{\Delta}}/\sqrt{d}<c{{\Delta}}/\sqrt{d},$ where ${{\Delta}}$ is the minimum separation between two centers. Once this density is available within small $L_{\infty}$ error, the local maxima can be approximately obtained using a robust randomized zeroth order concave optimization method developed in [3] started from all elements of a sufficiently rich $\ell_{\infty}^{d}$ net on $\bigcup_{i}B_{i}$ (which has polynomial size by Lemma 5), and the resulting centers are good approximations (i.e.  within $cd^{-\frac{5}{2}}$ of the true centers in Hausdorff distance) of the true centers with high probability. We then feed these $d^{-\frac{5}{2}}$ approximate centers into the iterative procedure developed by Regev and Vijayaraghavan in [13] and by Theorem 4.1 of the full version of that paper, a seed of this quality, suffices to produce in $\mathrm{poly}(k,d,{\mathbf{\delta}}^{-1})$ time, a set of centers that are within ${\mathbf{\delta}}$ of the true centers in Hausdorff distance.

The deconvolution is achieved by convolving the empirical measure $\mu_{e}$ with the Fourier transform of a certain carefully chosen $\hat{\zeta}$. The function $\hat{\zeta}$ is upto scalar multiplication, the reciprocal of a Gaussian with standard deviation $\sqrt{1-\bar{{{\Delta}}}^{2}}$ restricted to a ball of radius $(\sqrt{\log k}+\sqrt{d})\bar{{\Delta}}^{-1}$. It follows from Lemma 3, that the effect of this truncation on the deconvolution process can be controlled. The pointwise evaluation of the convolution is done using the Monte Carlo method. The truncation plays an important role, because without it, the reciprocal of a Gaussian would not be in $L_{p}$ for any $p\in[1,\infty]$, leading to difficulties in the Fourier analysis.

Suppose we are given independent, identically distributed samples $x_{1},\dots,x_{n}$ from a mixture $\mu$ of no more than $k$ of $d$-dimensional spherical gaussian distributions $\mu_{i}$ with variance $1$, such that the minimum $\ell_{2}$ distance between two distinct centers $y_{l}$ and $y_{j}$ is greater than $\sqrt{d}{{\Delta}}$ for some ${{\Delta}}\geq 2c_{0}$, where $c_{0}$ is certain universal constant in $(0,1)$. Suppose that $\mu:=\sum_{l=1}^{k_{0}}w_{l}\mu_{l}$, where $\mu_{l}$ are $d$-dimensional Gaussians with center at $y_{l}$, and $k_{0}\leq k$, and $(w_{1},\cdots,w_{k_{0}})$ is a known probability vector such that $w_{min}:=\min_{l\in[k_{0}]}w_{l}$ satisfies $w_{min}\geq\frac{c}{k}$.

Let $Y:=\{y_{1},\dots,y_{k_{0}}\}$ be the set of centers of the component gaussians in the mixture $\mu$:

Let us fix $1-\eta:=\frac{9}{10}$, to be the success probability we will require. This can be made $1-\eta_{0}$ that is arbitrarily close to $1$ by the following simple clustering technique in the metric space associated with Hausdorff distance, applied to the outputs of $100(\log\eta_{0}^{-1})$ runs of the algorithm.

Let $\mathcal{F}$ denote the fourier transform. For any function $f\in L^{2}(\mathbb{R}^{d})$, we have $\hat{f}=\mathcal{F}(f)$, where

By the fourier inversion formula,

where these are not necessarily lebesgue integrals, but need to be interpreted as improper integrals. Thus,

and

where the limit is taken in the $L^{2}$ sense.

Let

denote the standard gaussian in $d$ dimensions. Then, $\hat{\gamma}=\mathcal{F}(\gamma)=\gamma$. Let $\bar{{\Delta}}$ equal ${{\Delta}}/C_{3.2}$. Let

where $B_{2}(0,\bar{{\Delta}}^{-1})$ is the euclidean ball of radius $\bar{{\Delta}}^{-1}$ and $\mathbb{I}(S)$ the indicator function of $S$, and let $s=\mathcal{F}^{-1}(\hat{s})$. Let $\gamma_{\sigma}(w):=\sigma^{-d}\gamma(\sigma^{-1}w)$ denote the spherical gaussian whose one dimensional marginals have variance $\sigma^{2}$.

Let $\mu_{e}$ denote the uniform probability measure on $\{x_{1},\dots,x_{n}\}$. Let $\star$ denote convolution in $\mathbb{R}^{d}$. Note that the Fourier convolution identity is

Let $\hat{\zeta}:=\hat{\gamma}^{-1}\cdot\hat{s}$, and $\zeta=\mathcal{F}^{-1}(\hat{\zeta})$. We will recover the centers and weights of the gaussians from $\zeta\star\mu_{e}=(2\pi)^{\frac{d}{2}}\mathcal{F}^{-1}(\hat{\zeta}\cdot\hat{\mu}_{e})$. The heuristics are as follows.

Let $\nu$ denote the unique probability measure satisfying $\gamma\star\nu=\mu$. Thus

where ${\mathbf{\delta}}_{y_{j}}$ is a dirac delta supported on $y_{j}$. From this it follows, roughly speaking, that inside $\mbox{supp}(\hat{s})$, we get

pointwise, and this should (roughly) yield

On the other hand, notice that Lemma  3 shows that $s\star\nu\approx\gamma_{\bar{{\Delta}}}\star\nu$, and because the spikes of $\nu$ are approximately (up to scaling) the spikes of $\gamma_{\bar{\Delta}}\ast\nu$, we restrict to learning the spikes of $\gamma_{\bar{\Delta}}\ast\nu$ by accessing an approximation via $(2\pi)^{-\frac{d}{2}}\zeta\ast\mu_{e}$. For notational convenience, we write $\xi_{e}:=(2\pi)^{-\frac{d}{2}}\zeta\ast\mu_{e}$.

In order to formalize this, we proceed as follows. We will employ Hoeffding’s inequality for $\mathbb{C}$-valued random variables:

We write ${\mathcal{B}}:=\mbox{supp}(\hat{s})$, so that $z\in{\mathcal{B}}$ if and only if $z\in\mathbb{R}^{d}$ satisfies $||{\bar{\Delta}}z||\leq\sqrt{C_{3.5}\ln k}+\sqrt{d}$.

The following proposition allows us to construct a (random) black box oracle that outputs a good additive approximation of $\nu\star\gamma_{\bar{{\Delta}}}$ at any given point $x$.

We note that, for any $x\in\mathbb{R}^{d}$, one has

Writing $\xi(z):=\hat{s}(z)e^{\frac{||z||^{2}}{2}}\hat{\mu}_{e}(z)e^{ix\cdot z}$, one has