Ten Steps of EM Suffice for Mixtures of Two Gaussians

Introduction

The Expectation-Maximization (EM) algorithm [DLR77, Wu83, RW84] is one of the most widely used heuristics for maximizing likelihood in statistical models with latent variables. Consider a probability distribution $p_{\bm{\lambda}}$ sampling $(\bm{X},\bm{Z})$, where $\bm{X}$ is a vector of observable random variables, $\bm{Z}$ a vector of non-observable random variables and $\bm{\lambda}\in\bm{\Lambda}$ a vector of parameters. Given independent samples $\bm{x}_{1},\ldots,\bm{x}_{n}$ of the observed random variables, the goal of maximum likelihood estimation is to select $\bm{\lambda}\in\bm{\Lambda}$ maximizing the log-likelihood of the samples, namely $\sum_{i}\log p_{\bm{\lambda}}(\bm{x}_{i})$. Unfortunately, computing $p_{\bm{\lambda}}(\bm{x}_{i})$ involves summing $p_{\bm{\lambda}}(\bm{x}_{i},\bm{z}_{i})$ over all possible values of $\bm{z}_{i}$, which commonly results in a log-likelihood function that is non-convex with respect to $\bm{\lambda}$ and therefore hard to optimize. In this context, the EM algorithm proposes the following heuristic:

Start with an initial guess $\bm{\lambda}^{(0)}$ of the parameters.

(E-Step) For each sample $i$, compute the posterior $Q^{(t)}_{i}(\bm{z}):=p_{\bm{\lambda^{(t)}}}(\bm{Z}=\bm{z}|\bm{X}=\bm{x}_{i})$.

(M-Step) Set $\bm{\lambda^{(t+1)}}:=\arg\max_{\bm{\lambda}}\sum_{i}\sum_{\bm{z}}Q_{i}^{(t)}(\bm{z})\log{p_{\bm{\lambda}}(\bm{x}_{i},\bm{z})\over Q_{i}^{(t)}(\bm{z})}$.

Intuitively, the E-step of the algorithm uses the current guess of the parameters, $\bm{\lambda}^{(t)}$, to form beliefs, $Q^{(t)}_{i}$, about the state of the (non-observable) $\bm{Z}$ variables for each sample $i$. Then the M-step uses the new beliefs about the state of $\bm{Z}$ for each sample to maximize with respect to $\bm{\lambda}$ a lower bound on $\sum_{i}\log p_{\bm{\lambda}}(\bm{x}_{i})$. Indeed, by the concavity of the $\log$ function, the objective function used in the M-step of the algorithm is a lower bound on the true log-likelihood for all values of $\bm{\lambda}$, and it equals the true log-likelihood for $\bm{\lambda}=\bm{\lambda}^{(t)}$. From these observations, it follows that the above alternating procedure improves the true log-likelihood until convergence.

Despite its wide use and practical significance, little is known about whether and under what conditions EM converges to the true maximum likelihood estimator. A few works establish local convergence of the algorithm to stationary points of the log-likelihood function [Wu83, Tse04, CH08], and even fewer local convergence to the MLE [RW84, BWY17]. Besides local convergence, it is also known that badly initialized EM may settle far from the MLE both in parameter and in likelihood distance [Wu83]. The lack of theoretical understanding of the convergence properties of EM is intimately related to the non-convex nature of the optimization it performs.

Our paper aims to illuminate why EM works well in practice and develop techniques for understanding its behavior. We do so by analyzing one of the most basic and natural, yet still challenging, statistical models EM may be applied to, namely balanced mixtures of two multi-dimensional Gaussians with equal and known covariance matrices. In particular, we study the convergence of EM when applied to the following family of parametrized density functions:

where $\Sigma$ is a known covariance matrix, $(\bm{\mu}_{1},\bm{\mu}_{2})$ are unknown (vector) parameters, and $\mathcal{N}(\bm{\mu},\Sigma;\bm{x})$ represents the Gaussian density with mean $\bm{\mu}$ and covariance matrix $\Sigma$, i.e.

Our main contribution is to provide global convergence guarantees for EM applied to the above family of distributions. We establish our result for both the “population version” of the algorithm, and the finite-sample version, as described below.

To elucidate the optimization features of the algorithm and avoid analytical distractions arising due to sampling error, it has been standard practice in the literature of theoretical analyses of EM to consider the “population version” of the algorithm, where the EM iterations are performed assuming access to infinitely many samples from a distribution $p_{\bm{\mu}_{1},\bm{\mu}_{2}}$ as above. With infinitely many samples, we can identify the mean, ${\bm{\mu}_{1}+\bm{\mu}_{2}\over 2}$, of $p_{\bm{\mu}_{1},\bm{\mu}_{2}}$, and re-parametrize the density around the mean as follows:

We first study the convergence of EM when we perform iterations with respect to the parameter $\bm{\mu}$ of $p_{\bm{\mu}}(\bm{x})$ in (1.1). Starting with an initial guess $\bm{\lambda}^{(0)}$ for the unknown mean vector $\bm{\mu}$, the $t$-th iteration of EM amounts to the following update:

where we have compacted both the E- and M-step of EM into one update.

The intuition behind the EM update formula is as follows. First, we take expectations with respect to $\bm{x}\sim p_{\bm{\mu}}$ because we are studying the population version of EM, hence we assume access to infinitely many samples from $p_{\bm{\mu}}$. For each sample $\bm{x}$, the ratio $\frac{0.5\mathcal{N}(\bm{x};\bm{\lambda}^{(t)},\Sigma)}{p_{\bm{\lambda}^{(t)}}(\bm{x})}$ is our belief, at step $t$, that $\bm{x}$ was sampled from the first Gaussian component of $p_{\bm{\mu}}$, namely the one for which our current estimate of its mean vector is $\bm{\lambda}^{(t)}$. (The complementary probability is our present belief that $\bm{x}$ was sampled from the other Gaussian component.) Given these beliefs for all vectors $\bm{x}$, the update (1.2) is the result of the M-step of EM. Intuitively, our next guess $\bm{\lambda}^{(t+1)}$ for the mean vector of the first Gaussian component is a weighted combination over all samples $\bm{x}\sim p_{\bm{\mu}}$ where the weight of every $\bm{x}$ is our belief that it came from the first Gaussian component.

Our main result for population-EM is the following:

Whenever the initial guess $\bm{\lambda}^{(0)}$ is not equidistant to $\bm{\mu}$ and $-\bm{\mu}$, EM converges geometrically to either $\bm{\mu}$ or $-\bm{\mu}$, with convergence rate that improves as $t\rightarrow\infty$. We provide a simple, closed form expression of the convergence rate as a function of $\bm{\lambda}^{(t)}$ and $\bm{\mu}$. If the initial guess $\bm{\lambda}^{(0)}$ is equidistant to $\bm{\mu}$ and $-\bm{\mu}$, EM converges to the unstable fixed point $\bf 0$.

A formal statement is provided as Theorem 2 in Section 4. We start with the proof of the single-dimensional version, presented as Theorem 1 in Section 3. As a simple illustration of our result, we show in Section 5 that, in one dimension, when our original guess $\lambda^{(0)}=+\infty$ and the signal-to-noise ratio $\mu/\sigma=1$, $10$ steps of the EM algorithm result in $1\%$ error.

Despite the simplicity of the case we consider, no global convergence results were known prior to our work, even for the population EM. [BWY17] studied the same setting proving only local convergence, i.e. convergence only when the initial guess is close to the true parameters. They argue that the population EM update is contracting close to the true parameters. Unfortunately, the EM update is non-contracting outside a small neighborhood of the true parameters so this argument cannot be used for a global convergence guarantee.

In this work, we study the problem under arbitrary starting points and completely characterize the fixed points of EM. We show that other than a measure-zero subset of the space (namely points that are equidistant from the centers of the two Gaussians), any initialization of the EM algorithm converges to the true centers of the Gaussians, providing explicit bounds for the convergence rate. To achieve this, we follow an orthogonal approach to [BWY17]: Instead of trying to directly compute the number of steps required to reach convergence for a specific instance of the problem, we study the sensitivity of the EM iteration as the instance varies. The intuition is that if the EM update is sensitive to updating the instance, then changing the instance should also attract the update towards the changing instance; see Figure 1. We can use this, in turn, to argue that keeping the instance fixed, one EM update makes progress towards the true parameters. In particular, we gain a handle on the convergence rate of EM on all instances at once. This is quantified by Eq. (3.2).

Analysis of Finite-Sample EM for Mixtures of Two Gaussians.

Bootstrapping EM for Faster Convergence.

We note that, in multiple dimensions, care must be taken in initializing the EM algorithm, even in the infinite sample regime, as the convergence guarantee depends on the angle between the current iterate and the true mean vector. While a randomly chosen unit vector will have projection of $\Theta(1/\sqrt{d})$ in the direction of $\bm{\mu}$, we argue that we can boostrap EM to turn this projection larger than a constant. This allows us to work with similar convergence rates as in the single-dimensional case, namely only SNR (and not dimension) dependent. Our initialization procedure is described in Section 6.3.

EM can be boostrapped so that the number of iterations required to approximate $\bm{\mu}_{1}$ and $\bm{\mu}_{2}$ to within Mahalanobis distance $\epsilon$ depends logarithmically in the dimension.

Related Work on Learning Mixtures of Gaussians.

We have already outlined the literature on the Expectation-Maximization algorithm. Several results study its local convergence properties and there are known cases where badly initialized EM fails to converge. See above.

There is also a large body of literature on learning mixtures of Gaussians. A long line of work initiated by Dasgupta [Das99, AK01, VW04, AM05, KSV05, DS07, CR08, BV08, CDV09] provides rigorous guarantees on recovering the parameters of Gaussians in a mixture under separability assumptions, while later work [KMV10, MV10, BS10] has established guarantees under minimal information theoretic assumptions. More recent work [HP15] provides tight bounds on the number of samples necessary to recover the parameters of the Gaussians as well as improved algorithms, while another strand of the literature studies proper learning with improved running times and sample sizes [SOAJ14, DK14]. Finally, there has been work on methods exploiting general position assumptions or performing smoothed analysis [HK13, GHK15].

In practice, the most common algorithm for learning mixtures of Gaussians is the Expectation-Maximization algorithm, with the practical experience that it performs well in a broad range of scenarios despite the lack of theoretical guarantees. Recently, Balakrishnan, Wainwright and Yu [BWY17] studied the convergence of EM in the case of an equal-weight mixture of two Gaussians with the same and known covariance matrix, showing local convergence guarantees. In particular, they show that when EM is initialized close to the actual parameters, then it converges. In this work, we revisit the same setting considered by [BWY17] but establish global convergence guarantees. We show that, for any initialization of the parameters, the EM algorithm converges geometrically to the true parameters. We also provide a simple and explicit formula for the rate of convergence.

Concurrent and independent work by Xu, Hsu and Maleki [XHM16] has also provided global and geometric convergence guarantees for the same setting, as well as a slightly more general setting where the mean of the mixture is unknown, but they do not provide explicit convergence rates. They also do not provide an analysis of the finite-sample regime.

Preliminary Observations

In this section we illustrate some simple properties of the EM update (1.2) and simplify the formula. First, it is easy to see that plugging in the values $\bm{\lambda}\in\{-\bm{\mu},\bm{0},\bm{\mu}\}$ into $M(\bm{\lambda},\bm{\mu})$ results into

In particular, for all $\bm{\mu}$, these values are certainly fixed points of the EM iteration. Next, we rewrite $M(\bm{\lambda},\bm{\mu})$ as follows:

It is easy to observe that by symmetry this simplifies to

Simplifying common terms in the density functions $\mathcal{N}(\bm{x};\bm{\lambda},\Sigma)$, we get that

We thus get the following expression for the EM iteration

Single-dimensional Convergence

In the single dimensional case the EM algorithm takes the following form according to (2.2).

Observe that the function $M(\lambda,\mu)$ is increasing with respect to $\lambda$. Indeed the partial derivative of $M$ with respect to $\lambda$ is

which is strictly greater than zero since the $\tanh^{\prime}$ function is strictly positive.

We will show next that the fixed points we identified at (2.1) are the only fixed points of $M(\cdot,\mu)$. When initialized with $\lambda^{(0)}>0$ (resp. $\lambda^{(0)}<0$), the EM algorithm converges to $\mu>0$ (resp. to $-\mu<0$). The point $\lambda=0$ is an unstable fixed point.

In the single dimensional case, when $\lambda^{(0)},\mu>0$, the parameters $\lambda^{(t)}$ satisfy

Moreover $\kappa^{(t)}$ is a decreasing function of $t$.

For simplicity we will use $\lambda$ for $\lambda^{(t)}$, $\lambda^{\prime}$ for $\lambda^{(t+1)}$ and we will assume that $X\sim\mathcal{N}(0,\sigma^{2})$.

By a simple change of variables we can see that

The main idea is to use the Mean Value Theorem with respect to the second coordinate of the function $M$ on the interval $[\lambda,\mu]$.

But we know that $M(\lambda,\lambda)=\lambda$ and $M(\lambda,\mu)=\lambda^{\prime}$ and therefore we get

where we have used the fact that $\lambda^{\prime}<\mu$ which is comes from the fact that $M(\lambda,\mu)$ is increasing with respect to $\lambda$ and that $M(\mu,\mu)=\mu$.

The only thing that remains to complete our proof is to prove a lower bound of the partial derivative of $M$ with respect to $\mu$.

The first term is non-negative, Lemma 1. The second term is at least $1-\exp\left[-\frac{\min(\xi,\lambda)\cdot\xi}{2\sigma^{2}}\right]$, Lemma 2 and the theorem follows. ∎

But now we can see that since $\tanh^{\prime}$ is an even function and since for any $y>0$ we have $\exp{\left(-\frac{(y-\alpha)^{2}}{2\sigma^{2}}\right)}\geq\exp{\left(-\frac{(-y-\alpha)^{2}}{2\sigma^{2}}\right)}$ then

Multi-dimensional Convergence

In the multidimensional case, the EM algorithm takes the form of (2.2). In this case, we will quantify our approximation guarantees using the Mahalanobis distance $\left\|\cdot\right\|_{\Sigma}$ between vectors with respect to matrix $\Sigma$, defined as follows:

We will show that the fixed points identified in (2.1) are the only fixed points of $M(\cdot,\bm{\mu})$. When initialized with $\bm{\lambda}^{(0)}$ such that $\left\|\bm{\lambda}^{(0)}-\bm{\mu}\right\|_{\Sigma}<\left\|\bm{\lambda}^{(0)}+\bm{\mu}\right\|_{\Sigma}$ (resp. $\left\|\bm{\lambda}^{(0)}-\bm{\mu}\right\|_{\Sigma}>\left\|\bm{\lambda}^{(0)}+\bm{\mu}\right\|_{\Sigma}$), the EM algorithm converges to $\bm{\mu}$ (resp. to $-\bm{\mu}$). The algorithm converges to $\bm{\lambda}=\bm{0}$ when initialized with $\left\|\bm{\lambda}^{(0)}-\bm{\mu}\right\|_{\Sigma}=\left\|\bm{\lambda}^{(0)}+\bm{\mu}\right\|_{\Sigma}$. In particular,

Whenever $\left\|\bm{\lambda}^{(0)}-\bm{\mu}\right\|_{\Sigma}<\left\|\bm{\lambda}^{(0)}+\bm{\mu}\right\|_{\Sigma}$, i.e. the initial guess is closer to $\bm{\mu}$ than $-\bm{\mu}$, the estimates $\bm{\lambda}^{(t)}$ of the EM algorithm satisfy

Moreover, $\kappa^{(t)}$ is a decreasing function of $t$. The symmetric things hold when $\left\|\bm{\lambda}^{(0)}-\bm{\mu}\right\|_{\Sigma}>\left\|\bm{\lambda}^{(0)}+\bm{\mu}\right\|_{\Sigma}$. When the initial guess is equidistant to $\bm{\mu}$ and $-\bm{\mu}$, then $\bm{\lambda}^{(t)}=\bm{0}$ for all $t>0$.

For simplicity we will use $\bm{\lambda}$ for $\bm{\lambda}^{(t)}$, $\bm{\lambda}^{\prime}$ for $\bm{\lambda}^{(t+1)}$.

By applying the following change of variables $\bm{\lambda}\leftarrow\Sigma^{-1/2}\bm{\lambda}$ and $\bm{\mu}\leftarrow\Sigma^{-1/2}\bm{\mu}$ we may assume that $\Sigma=I$ where $I$ is the identity matrix. Therefore the iteration of EM becomes

Since the Normal distribution is rotation invariant we can equivalently write:

We now consider different cases for $i$ to further simplify Equation (4.2).

Let $\kappa$ as defined before and using Lemma 2 we get that

We can now bound the distance of $\bm{\lambda}^{\prime}$ from $\bm{\mu}$:

We now have to prove that this convergence rate $\kappa$ decreases as the iterations increase. This is implied by the following lemmas which show that $\min(\langle\hat{\bm{\lambda}},\bm{\lambda}\rangle,\langle\hat{\bm{\lambda}},\bm{\mu}\rangle)\leq\min(\langle\hat{\bm{\lambda}}^{\prime},\bm{\lambda}^{\prime}\rangle,\langle\hat{\bm{\lambda}}^{\prime},\bm{\mu}\rangle)$

If $\left\|\bm{\lambda}\right\|\geq\langle\hat{\bm{\lambda}},\bm{\mu}\rangle$ then $\langle\hat{\bm{\lambda}},\bm{\mu}\rangle\leq\left\|\bm{\lambda}^{\prime}\right\|$ and $\langle\hat{\bm{\lambda}},\bm{\mu}\rangle\leq\langle\hat{\bm{\lambda}^{\prime}},\bm{\mu}\rangle$.

The analysis above implies that $\bm{\lambda}^{\prime}$ can be written in the form $\bm{\lambda}^{\prime}=\alpha\cdot\hat{\bm{\lambda}}+\beta\cdot\hat{\bm{\lambda}}^{\bot}$, where $\langle\hat{\bm{\lambda}},\bm{\mu}\rangle\leq\alpha\leq\left\|\bm{\lambda}\right\|$ and $0\leq\beta\leq\langle\hat{\bm{\lambda}}^{\bot},\bm{\mu}\rangle$. It is easy to see that the first inequality holds since $\left\|\bm{\lambda}^{\prime}\right\|\geq\alpha\geq\langle\hat{\bm{\lambda}},\bm{\mu}\rangle$. For the second, we write $\langle\hat{\bm{\lambda}^{\prime}},\bm{\mu}\rangle$ as:

where we used the fact that $\frac{\langle\hat{\bm{\lambda}}^{\bot},\bm{\mu}\rangle}{\langle\hat{\bm{\lambda}},\bm{\mu}\rangle}\geq\frac{\beta}{\alpha}$ which follows by the bounds on $\alpha$ and $\beta$. ∎

If $\left\|\bm{\lambda}\right\|\leq\langle\hat{\bm{\lambda}},\bm{\mu}\rangle$ then $\left\|\bm{\lambda}\right\|\leq\left\|\bm{\lambda}^{\prime}\right\|\leq\langle\hat{\bm{\lambda}^{\prime}},\bm{\mu}\rangle$.

We have that $\bm{\lambda}^{\prime}=\alpha\cdot\hat{\bm{\lambda}}+\beta\cdot\hat{\bm{\lambda}}^{\bot}$, where $\left\|\bm{\lambda}\right\|\leq\alpha\leq\langle\hat{\bm{\lambda}},\bm{\mu}\rangle$ and $0\leq\beta\leq\langle\hat{\bm{\lambda}}^{\bot},\bm{\mu}\rangle$. We also have $\langle\bm{\lambda}^{\prime},\bm{\mu}\rangle=\alpha\langle\hat{\bm{\lambda}},\bm{\mu}\rangle+\beta\langle\hat{\bm{\lambda}}^{\bot},\bm{\mu}\rangle\geq\alpha^{2}+\beta^{2}=\left\|\bm{\lambda}^{\prime}\right\|^{2}\geq\alpha^{2}\geq\left\|\bm{\lambda}\right\|^{2}$ so the lemma follows. ∎

Finally substituting back in the basis that we started before changing coordinates to make the covariance matrix identity we get the result as stated at the theorem. ∎

An Illustration of the Speed of Convergence

Using our results in the previous sections we can calculate explicit speeds of convergence of EM to its fixed points. In this section, we present some results with this flavor. For simplicity, we start with single dimensional case, and discuss the multi-dimensional case in the end of this section.

Let us consider a mixture of two single-dimensional Gaussians whose signal-to-noise ratio $\eta=\mu/\sigma$ is equal to $1$. There is nothing special about the value of $1$, except that it is a difficult case to consider since the Gaussian components are not separated, as shown in Figure 3.

When the SNR is larger, the numbers presented below still hold and in reality the convergence is even faster. When the SNR is even smaller than one, the numbers change, but gracefully, and they can be calculated in a similar fashion.

We will also assume a completely agnostic initialization of EM, setting $\lambda^{(0)}\rightarrow+\infty$.In the multi-dimensional setting, this would corrrespond to a very large magnitude $\bm{\lambda}^{(0)}$ chosen in a random direction. To analyze the speed of convergence of EM to its fixed point $\mu$, we first make the observation that in one step we already get to $\lambda^{(1)}\leq\mu+\sigma\sqrt{2\over\pi}$. To see this we can plug $\lambda^{(0)}\rightarrow\infty$ into equation (3.1) to get:

which equals the mean of the Folded Normal Distribution. A well-known bound for this mean is $\mu+\sqrt{2\over\pi}\sigma$. Therefore the distance from the true mean after one step is $\left|\lambda^{(1)}-\mu\right|\leq\sqrt{2\over\pi}\sigma$.

Now, using Theorem 1, we conclude that in all subsequent steps the distance to $\mu$ shrinks by a factor of at least $e^{+1/2}$. This means that, if we want to estimate $\mu$ to within additive error $1\%\sigma$, then we need to run EM for at most $\lceil 2\cdot\ln 100+\ln{2\over\pi}\rceil=9$ additional steps. Accounting for the first step, $10$ iterations of the EM algorithm in total suffice to get to within error $1\%$, even when our initial guess of the mean is infinitely away from the true value!

The analysis given above for $SNR=1$, generalizes to arbitrary $SNR$s, at the cost of a factor of $O(1/SNR^{2})$ in the number of iterations. It also generalizes to obtain an arbitrary approximation $\epsilon$, at a cost of a factor of $O(\log 1/\epsilon)$. In multiple dimensions, we could run EM from a random initialization. The number of iterations for an approximation of $\epsilon$ in Mahalanobis distance would depend on the angle of the initial iterate with $\bm{\mu}$. Under a random initialization, the cosine of that angle is expected to be $\Theta(1/\sqrt{d})$, resulting in a total $O(d/SNR^{2}\cdot\log(1/\epsilon))$ number of iterations. We show that we can boostrap EM to obtain a better initialization, starting from a random one, improving the angle to $\Omega(1)$, after $O(\log d/SNR^{2})$ iterations. With a constant angle, EM takes $O(1/SNR^{2}\cdot\log(1/\epsilon))$ iterations to give $\epsilon$ error, as in the single dimension, overall improving exponentially the dependence on $d$. We describe our bootstrapping operation in the context of our analysis of the finite sample EM in Section 6.3.

Sample Based Model

The main goal of this section is to prove convergence guarantees for the EM algorithm, when we have a finite sample. Similarly to the Section 4 we willquantify our approximation guarantees using the Mahalanobis distance $\left\|\cdot\right\|_{\Sigma}$ between two vectors with respect to matrix $\Sigma$, which we remind is defined as follows:

Also for the simplicity of the notation, it is useful to define the Mahalanobis inner product between two vectors with respect to matric $\Sigma$ as follows:

Towards our goal, we encounter two challenges.

The first is that we cannot assume that we exactly know the mean of the mixture distribution $(\bm{\mu}_{1}+\bm{\mu}_{2})/2$. Our only access to this mean is via samples. We therefore use samples to estimate it. Then, we translate the origin to our estimate, and write the EM iteration for finding the mean of one of the two mixture components with respect to this origin. Given the error incurred in the approximation of $(\bm{\mu}_{1}+\bm{\mu}_{2})/2$, we propose to stabilize the sample-based EM iteration by including in the sample for each sampled point $\bm{x}_{i}$ its symmetric point $-\bm{x}_{i}$. This is the sample based version that we analyze, although our analysis goes through without this stabilization.

We use $\bm{\delta}$ to refer to the error in our estimation $\bm{c}-\frac{\bm{\mu}_{1}+\bm{\mu}_{2}}{2}$.

Combining the above steps all together we get our main theorem for this section.

We start now proving the lemmas for each of the steps described above.

We start by making the transformation $x\mapsto\Sigma^{-1/2}x$ to the space, so that the covariance matrix is the identity $I$. Finally we will take back this transformation and the Euclidean norm becomes the corresponding Mahalanobis.

We will get an estimate $\bm{c}$ of the mean $(\bm{\mu}_{1}+\bm{\mu}_{2})/2$ by drawing $O(d/\varepsilon^{2})$ samples from the mixture and working in each axis direction separately. We will compute the estimate $\hat{c}_{i}$ in axis direction $\bm{e}_{i}$ as the average of the first and third quartile of the empirical distribution given by the samples. It suffices to show that with high probability every quartile is at most $\frac{\varepsilon}{\sqrt{d}}$ away from the true quartile of the distribution $p_{\bm{\mu}_{1},\bm{\mu}_{2}}\cdot\bm{e}_{i}$.

Let $p_{-\mu,\mu}$ be the mixture of Gaussian distributions obtained by centering $p_{\bm{\mu}_{1},\bm{\mu}_{2}}\cdot\bm{e}_{i}$ around .

The cumulative distribution $F$ of $p_{-\mu,\mu}$ is given by $F(x)=\frac{1}{2}\Phi\left(x+\mu\right)+\frac{1}{2}\Phi\left(x-\mu\right)$. By the DKW inequality [DL12] with $n$ samples we have that the empirical distribution $F_{n}$ satisfies:

In particular, for $x$ such that $F_{n}(x)=1/4$, with high probability $F(x)=1/4\pm\varepsilon/\sqrt{d}$. Moreover, for $x^{*}$ such that $F(x^{*})=1/4$, we have that $|x-x^{*}|\leq\frac{\varepsilon/\sqrt{d}}{\min_{\xi\in[x,x^{*}]}F^{\prime}(\xi)}$ by the mean value theorem.

Since for $\xi\in[x,x^{*}]$, $F(\xi)\in[1/4-\varepsilon/\sqrt{d},1/4+\varepsilon/\sqrt{d}]$, this implies that $\Phi\left(\xi+\mu\right)\in[1/4-\varepsilon/\sqrt{d},1/2+2\varepsilon/\sqrt{d}]$ and thus $F^{\prime}(\xi)\geq\frac{1}{2}\Phi^{\prime}\left(\xi+\mu\right)$. But $\Phi^{\prime}(z)\geq\frac{1}{5}$ when $\Phi(z)\in[1/5,4/5]$ which shows that $|x-x^{*}|=O(\varepsilon/\sqrt{d})$. Similarly, this holds for the 3rd quartile as well and thus the same bound holds for the mean as well of the two quartiles with probability $1-2e^{{-2n\varepsilon^{2}/d}}$. Setting $n=O(d\log d/\varepsilon^{2})$, we get that the bound is violated with probability $O(1/d)$. Taking a union bound for all $d$ axis directions, we get that for all $i$, $|\bm{\delta}_{i}|=O(\varepsilon/\sqrt{d})$ with constant probability. This implies that $\left\|\bm{\delta}\right\|\leq\varepsilon$ with constant probability. By using a factor of $\log 1/\eta$ more samples it is easy to see that the error probability reduces to $\eta$. ∎

2 Sample Based EM Iteration

Using the estimation of the center $\bm{c}$ that we found in the previous section, we translate all our data and parameters so that $\bm{c}\mapsto\bm{0}$. After this centering, the parameters $\bm{\mu}_{1}$, $\bm{\mu}_{2}$ become $\bm{\mu}+\bm{\delta}$, $-\bm{\mu}+\bm{\delta}$ and the covariance matrix remains the same.

3 Initialization of EM

We rewrite the sample based EM iteration is the following form

where $\hat{\Sigma}_{p}$ is the empirical covariance matrix of the mixture $p_{\bm{\mu}+\bm{\delta},-\bm{\mu}+\bm{\delta}}$. Now the intuition suggests that the direction of the maximum eigenvector of the matrix $\hat{\Sigma}_{p}\Sigma^{-1}$ is a direction that is spanned by $\bm{\mu}+\bm{\delta}$, $-\bm{\mu}+\bm{\delta}$ and hence a direction with small angle to the direction of $\bm{\mu}$. Also observe that this approximate EM iteration is actually an iteration of the power method for the matrix $\hat{\Sigma}_{p}\Sigma^{-1}$! Because of the efficiency of the power method, we expect that after a few steps this iteration will find a direction $\hat{\bm{\lambda}}$ such that the inner product $\langle\hat{\bm{\mu}},\hat{\bm{\lambda}}\rangle$ is large enough. This direction is a good initialization for the EM algorithm as we will see in the next section. We now formaly demonstrate the intuition we described for the approximate EM update.

We start by bounding the error we introduce by replacing the $\tanh$ function with its linear approximation. It is very easy to see that $\left|\tanh^{\prime\prime}(x)\right|\leq 1$. Now by Taylor expansion of $\tanh$ around we get that

At this point the calculations become much easier if we assume that we have already done the mapping $x\mapsto\Sigma^{-1/2}x$ and when we are done we will take the inverse mapping and get the result in Mahalanobis distance. Equation (6.4) suggests that it suffices to analyze the convergence of the power method given by the following equation.

where $\hat{\Sigma}_{p}$ is the empirical covariance matrix of the mixture $p_{\bm{\mu}+\bm{\delta},-\bm{\mu}+\bm{\delta}}$. Before analyzing (6.5) lets see what happens if instead of the empirical covariance $\hat{\Sigma}_{p}$ we had the actual covariance $\Sigma_{p}$ of the mixture distribution $p_{\bm{\mu}+\bm{\delta},-\bm{\mu}+\bm{\delta}}$. Then the iteration would be $\bm{\tau}^{(t+1)}=\Sigma_{p}\bm{\tau}^{(t)}$. The covariance matrix of the mixture is

Therefore the principal eigenvector of $\Sigma_{p}$ is $\hat{\bm{\mu}}=\frac{\bm{\mu}}{\left\|\bm{\mu}\right\|}$ with eigenvalue $1+\mu^{2}$, where $\mu=\left\|\bm{\mu}\right\|$. All the other eigenvectors have eigenvalue $1$ and therefore the ratio of the largest to the lowest eigenvalue is $\rho=1+\mu^{2}$.

To get the corresponding properties of $\hat{\Sigma}_{p}$ we observe that each $\bm{x}_{i}$ can be written as

where $\bm{y}_{i}$ is distributed as $\mathcal{N}(0,I)$ and $z_{i}$ is a Rademacher indicator variable that shows whether $\bm{x}_{i}$ is coming from the distribution $\mathcal{N}(\bm{\mu}+\bm{\delta},I)$ or $\mathcal{N}(\bm{\mu}-\bm{\delta},I)$. We notice that $\bm{y}_{i}$ and $z_{i}$ are independent. We can now rewrite $\hat{\Sigma}_{p}$ in terms of $\bm{y}_{i}$ and $z_{i}$.

We describe now a sketch of the proof of Lemma 6, using the following Lemma 25 from [DKT15].

For all $i\in\{1,\dots d\}$, $\left|\left(\frac{\bm{v}_{i}}{\sqrt{\lambda_{i}}}\right)^{T}\left(\Sigma-\hat{\Sigma}\right)\left(\frac{\bm{v}_{i}}{\sqrt{\lambda_{i}}}\right)\right|\leq\epsilon$,

For all $i,j\in\{1,\dots d\}$, $\left|\left(\frac{\bm{v}_{i}}{\sqrt{\lambda_{i}}}+\frac{\bm{v}_{j}}{\sqrt{\lambda_{j}}}\right)^{T}\left(\Sigma-\hat{\Sigma}\right)\left(\frac{\bm{v}_{i}}{\sqrt{\lambda_{i}}}+\frac{\bm{v}_{j}}{\sqrt{\lambda_{j}}}\right)\right|\leq 4\epsilon$.

Now we consider the quantity $\bm{e}_{j}^{T}\left(\frac{1}{n}\sum_{i=1}^{n}\bm{y}_{i}\right)$, where $\bm{e}_{j}$ is the unit $j$th vector. This is equivalent with having $\frac{1}{n}\sum_{i=1}^{n}y_{i}$ where $y_{i}$ is drawn from $\mathcal{N}(0,1)$. So we have that

Now by doing a union bound over all $\bm{e}_{j}$ we get that

For $\Sigma_{2}$, it is easy to observe that in $\hat{\Sigma}_{2}=\frac{1}{n}\sum_{i=1}^{n}z_{i}\left(\bm{y}_{i}\bm{\delta}^{T}+\bm{\delta}\bm{y}_{i}^{T}\right)$, $z_{i}$ are independent from $\bm{y}_{i}$ and so the product $z_{i}\bm{y}_{i}$ is a sample from standard multinormal distribution $\mathcal{N}(0,I)$. Therefore we can substitute $z_{i}\bm{y}_{i}$ with just $\bm{y}_{i}$. Now using exactly the same analysis as in Lemma 7 we can prove the following lemma.

Now we calculate the variance of the samples in the direction of $\bm{\mu}$. We have

The principal eigenvector $\hat{\bm{v}}$ of $\hat{\Sigma}_{p}$ has to satisfy

Now using an $\varepsilon$ such that $\varepsilon\leq\mu/10$ we have that the above implies $a\geq 3/4$. This proves the following Proposition that we use in the analysis of the initialization step.

The last thing we need to prove to complete the analysis of the initialization is the gap between the first and the second eigenvalue of $\hat{\Sigma}_{p}$. We want to use this gap for the analysis of the convergence rate of the power method iteration that takes place in the first steps. Given Proposition refprop:principalVectorSp, we have that any direction $\hat{v}$ except from the principal one, has $\langle\hat{\bm{v}},\hat{\bm{\mu}}\rangle\leq\frac{1}{4}$ with high probability. Let $\hat{\bm{v}}^{\prime}$ be the eigenvector that corresponds to the second maximum eigenvalue of $\hat{\Sigma}_{p}$. Let $a=\hat{\bm{v}}^{\prime T}\bm{\mu}$, using Lemmas 6, 8, 7 and (6.7) we have that

On the other hand based on the fact that $\hat{\bm{v}}^{T}\hat{\Sigma}_{p}\hat{\bm{v}}\geq\hat{\bm{\mu}}^{T}\hat{\Sigma}_{p}\hat{\bm{\mu}}$ we conclude that

Using also the hypothesis that $0\leq\varepsilon\leq\mu/10$ we get that

It is well known and easy to prove that if we choose a random vector $\bm{\lambda}^{(0)}$ uniformly from the half unit sphere, defined by $\hat{\bm{\mu}}$, then we will have that $\langle\bm{\lambda}^{(0)},\hat{\bm{\mu}}\rangle\geq 1/\sqrt{d}$ with high probability.

But $k+1$ is $O(\log d/\log\hat{\rho})$ and therefore

Finally since $\hat{\rho}$ is $\min\left(1+\frac{\mu^{2}}{4},2\right)$ and also $\varepsilon\leq\mu/10$ we only need to set $\epsilon$ polynomially with respect to $\varepsilon^{2}/d$ and we will get that after the first $O(\log d/\min(\mu^{2},1))$ iterations it holds that

for some $\eta\leq\varepsilon/6$. Which implies that

This proves the following lemma and completes the proof of the initialization.

4 Finite-Sample EM Analysis

We initialize EM at $\bm{\lambda}^{(0)}=M\bm{\lambda}$, where $\bm{\lambda}$ is the point from Lemma 9, for some large constant $M$.It is easy to find such constant by getting a small number of samples and keeping the one that has maximum magnitude. With this initialization, we run EM for $t=O((1/\mu^{2})\log(d/\varepsilon))$ steps using $O((d/\varepsilon^{2}\mu^{4})\log 1/\delta)$ samples at each step, where for ease of notation we have set $\mu=SNR\equiv\left\|\bm{\mu}\right\|_{\Sigma}$.

To study our sample-based EM iteration (6.1) we will relate its progress to an appropriate population EM iteration. Note that this iteration differs from the population EM iteration that we discussed in Section 2 and analyzed in Section 4. The reason is that we have incurred an error $\bm{\delta}$ in the estimation of the mean of the distribution in Section 6.1. With respect to our estimated mean centering, the true means of the two Gaussian components are $\bm{\mu}+\bm{\delta}$, $-\bm{\mu}+\bm{\delta}$ rather than $\bm{\mu}$ and $-\bm{\mu}$. Another source of discrepancy comes from the fact that we included for each point $\bm{x}_{i}$ in our sample its symmetric point $-\bm{x}_{i}$. This implies that each $\bm{x}_{i}$ is coming with probability $1/2$ from the mixture $p_{\bm{\mu}+\bm{\delta},-\bm{\mu}-\bm{\delta}}$ and with probability $1/2$ from the mixture $p_{\bm{\mu}-\bm{\delta},-\bm{\mu}+\bm{\delta}}$. Given this, using again the same operations as in Section 2, we have that the corresponding population iteration, denoted by $\bm{\lambda}^{(t)}$, is

Our proof follows two steps illustrated in Figures 2 and LABEL:fig:graph3:

Step 1: First, we relate the population EM iteration defined by (6.8) to the vanilla population EM iteration defined by (2.2);

Step 2: Then, we related the population EM iteration defined by (6.8) to the sample-based iteration.

To analyze the convergence of (6.8), we use Theorem 2 for every component of the mixture. More precisely, let $\bm{\lambda}_{1}^{(t)}$ and $\bm{\lambda}_{2}^{(t)}$ be

We are ready now to bound the convergence of $\bm{\lambda}^{(t)}=(\bm{\lambda}_{1}^{(t)}+\bm{\lambda}_{2}^{(t)})/2$

We can use the analysis of Section 4 to see that $\kappa^{\prime(t)}\leq\kappa^{\prime(0)}$. Therefore we also know that the population iteration satisfies:

where for simplicity we let $\kappa=\kappa^{\prime(0)}$. Now we have $\bm{\lambda}^{(0)}=M\bm{\lambda}$. Also $M$ is larger than $\mu$ and because of Lemma 9 $\langle\hat{\bm{\lambda}},\hat{\bm{\mu}}\rangle_{\Sigma}\geq 1/2$ with high probability, where $\mu=\left\|\bm{\mu}\right\|_{\Sigma}$. Therefore

and since $\varepsilon\leq\mu/10$ we have that

Also we have that $\left\|\bm{\delta}\right\|_{\Sigma}\leq\varepsilon$ and therefore (6.11) becomes

Step 2:

Now we define $\mu_{1}=\langle\bm{\mu},\bm{v}_{1}\rangle$, $\delta_{1}=\langle\bm{\delta},\bm{v}_{1}\rangle$ and we have

where $x_{i}$ is distributed as $p_{\mu_{1}+\delta_{1},\mu_{1}-\delta_{1}}$ and $\lambda=\left\|\bm{\lambda}\right\|$. For simplicity we refer to $p_{\mu_{1}+\delta_{1},\mu_{1}-\delta_{1}}$ as $\mathcal{D}_{1}$. Our goal is to bound the following probability

to do so we use the general large deviation technique. Because of symmetry of $\mathcal{D}_{1}$, we have that the above probability is equal twice the probability

Now we use the following Lemma 2.1 of [Wai15].

where $\bm{x},\bm{y}\sim\mathcal{N}(0,I)$ are standard multivariate Gaussian, and independent.

Now for the second term of (6.17) we notice that as we present in Section 3

which imlies using the mean value theorem that

Putting all together to (6.17) we have that

As before we define $\mu_{2}=\langle\bm{\mu},\bm{v}_{2}\rangle$, $\delta_{2}=\langle\bm{\delta},\bm{v}_{2}\rangle$ and we have

where $x_{i}$ is distributed as $p_{\mu_{1}+\delta_{1},\mu_{1}-\delta_{1}}$, $y_{i}$ is distributed as $p_{\mu_{2}+\delta_{2},\mu_{2}-\delta_{2}}$ and $\lambda=\left\|\bm{\lambda}\right\|$. For simplicity we refer to $p_{\mu_{2}+\delta_{2},\mu_{2}-\delta_{2}}$ as $\mathcal{D}_{2}$. Our goal is to bound the following probability

to do so we use the general large deviation technique. Using the symmetry of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ we have that the above probability is equal twice the

Using the fact that because of the initialization of EM $\lambda\geq\mu_{1}$ and by assumption $\delta_{1}\leq\mu_{1}$ and also let $\alpha=\min(1,\mu_{1})$ we have to bound the quantity

Now using a simple Taylor expansion used in the proof of the Hoeffding bound we get that the first term of (6.22) is less than or equal to

Now for the second term of (6.22) we have that

Putting all together to (6.22) we have that

For any $i\geq 3$ we follow the same analysis as for the bound (6.23) but because of the definition of the basis $\{\bm{v}_{1},\dots,\bm{v}_{d}\}$ we have that $\mu_{i}=\langle\bm{\mu},\bm{v}_{i}\rangle=0$ and therefore

Finally if we combine (6.18), (6.23) and (6.24) using the observation that $\left\|\bm{\lambda}^{(t+1)}-\bm{\mu}\right\|\geq\left|\mu_{2}\right|$ we get that

Proof of Theorem 3:

Acknowledgements

We thank Sham Kakade for suggesting the problem to us, and for initial discussions. The authors were supported by NSF Awards CCF-0953960 (CAREER), CCF-1551875, CCF-1617730, and CCF-1650733, ONR Grant N00014-12-1-0999, and a Microsoft Faculty Fellowship.

References