Comparing EM with GD in Mixture Models of Two Components

TL;DR

This paper compares the behavior of EM and gradient descent in two-component mixture models, showing EM generally escapes problematic regions faster and is less prone to trapping in bad local minima than GD.

Contribution

It provides a detailed analysis of how EM and GD differ in escaping one-cluster regions in Gaussian and Bernoulli mixture models, highlighting EM's advantages.

Findings

EM escapes one-cluster regions exponentially fast in Gaussian mixtures

GD escapes linearly fast and can get trapped in stable regions

EM can escape certain bad local minima that trap GD

Abstract

The expectation-maximization (EM) algorithm has been widely used in minimizing the negative log likelihood (also known as cross entropy) of mixture models. However, little is understood about the goodness of the fixed points it converges to. In this paper, we study the regions where one component is missing in two-component mixture models, which we call one-cluster regions. We analyze the propensity of such regions to trap EM and gradient descent (GD) for mixtures of two Gaussians and mixtures of two Bernoullis. In the case of Gaussian mixtures, EM escapes one-cluster regions exponentially fast, while GD escapes them linearly fast. In the case of mixtures of Bernoullis, we find that there exist one-cluster regions that are stable for GD and therefore trap GD, but those regions are unstable for EM, allowing EM to escape. Those regions are local minima that appear universally in experiments and can be arbitrarily bad. This work implies that EM is less likely than GD to converge to certain bad local optima in mixture models.

Figures

Tables

Table 1. Table 1: Average fraction of times (as a percentage) when GD converges to a k 𝑘 k -cluster point with k < m 𝑘 𝑚 k<m from random initialization. Infinite sample case.

	$D=2$	$D=3$	$D=4$	$D=5$	$D=6$
$m=2$	21.7	13.3	10.0	21.7	18.3
$m=3$	46.7	15.0	56.7	26.7	26.7
$m=4$	18.3	26.7	30.0	36.7	48.3
$m=5$	48.3	56.7	50.0	58.3	63.3
$m=6$	35.0	46.7	68.3	58.3	53.3

Table 2. Table 2: Average fraction of times (as a percentage) when GD converges to a k 𝑘 k -cluster point with k < m 𝑘 𝑚 k<m from random initialization. Finite sample case.

	$D=2$	$D=3$	$D=4$	$D=5$	$D=6$
$m=2$	10.0	5.0	0.0	13.3	10.0
$m=3$	10.0	13.3	21.6	36.7	40.0
$m=4$	33.3	35.0	10.0	36.7	56.7
$m=5$	33.3	46.7	25.0	40.0	70.0
$m=6$	43.3	50.0	58.3	60.0	80.0

Table 3. Table 3: Average ratio of the likelihood (as a percentage) of the k 𝑘 k -cluster points ( k < m 𝑘 𝑚 k<m ) that GD converges to, to the optimal likelihood. Infinite sample case.

	$D=4$	$D=5$	$D=6$	$D=8$	$D=10$
$m=2$	95.9	89.2	91.5	75.2	72.0
$m=3$	98.7	97.9	97.1	91.9	83.4
$m=4$	99.5	98.1	99.1	95.1	88.4
$m=5$	99.9	99.8	98.9	94.9	90.3
$m=6$	100.0	99.7	99.6	95.9	92.6

Table 4. Table 4: Worst-case ratio of the likelihood (as a percentage) of the k 𝑘 k -cluster points ( k < m 𝑘 𝑚 k<m ) that GD converges to, to the optimal likelihood. Infinite sample case.

	$D=4$	$D=5$	$D=6$	$D=8$	$D=10$
$m=2$	83.3	74.3	67.7	59.1	52.7
$m=3$	91.3	91.3	82.6	73.2	60.6
$m=4$	90.8	83.1	92.0	83.3	71.5
$m=5$	96.2	94.4	94.5	80.1	82.3
$m=6$	98.6	98.4	97.1	94.5	70.3

Equations

$$p(\mathbf{x})=\sum _{c=1}^m{\pi }_{c}f(\mathbf{x}|{\mathbf{\mu }}_c).$$

$$\ell =-\text{mathds}E[\log p(\mathbf{x})],$$

$$f(\mathbf{x}|{\mathbf{\mu }}_c)=\mathcal{N}(\mathbf{x}|{\mathbf{\mu }}_c,I),$$

$$f(\mathbf{x}|{\mathbf{\mu }}_c)=\prod _{i=1}^{D}B(x_i|{\mu }_{c,i})=B(\mathbf{x}|{\mathbf{\mu }}_c),$$

$$\text{mathds}{E}_{\mathbf{x}\sim B(\mathbf{x}|{\mathbf{\mu }}_c)}[\mathbf{x}]={\mathbf{\mu }}_c,$$

$${\gamma }_c(\mathbf{x}):=\frac{f(\mathbf{x}|{\mathbf{\mu }}_c)}{p(\mathbf{x})}.$$

$$M({\mathbf{\mu }}_c):=\frac{\text{mathds}E[\mathbf{x}{\gamma }_c(\mathbf{x})]}{\text{mathds}E[{\gamma }_c(\mathbf{x})]},M({\pi }_c):={\pi }_c\text{mathds}E[{\gamma }_c(\mathbf{x})].$$

$$Z_c=\int {\tilde{q}}_c(\mathbf{x})d\mathbf{x},q_c(\mathbf{x})=\frac{{\tilde{q}}_c(\mathbf{x})}{Z_c},$$

$$M({\mathbf{\mu }}_c):=\text{mathds}{E}_{\mathbf{x}\sim q_c(\mathbf{x})}[\mathbf{x}],M({\pi }_c):={\pi }_c{Z}_c.$$

$$\frac{\partial \ell }{\partial {\pi }_c}=-\text{mathds}E[{\gamma }_c(\mathbf{x})]=-Z_c.$$

$$\frac{\partial \ell }{\partial {\mathbf{\mu }}_c}$$

$$=$$

$$\frac{\partial \ell }{\partial {\mu }_{cj}}=-{\pi }_c\text{mathds}E\left[{\gamma }_c(\mathbf{x})\frac{(-1{)}^{1-x_j}}{B(x_j|{\mu }_{cj})}\right].$$

$$P_{{\mathbf{\mu }}_c}({\mathbf{\mu }}_c)=\min \{\max \{{\mathbf{\mu }}_c,0\},1\},$$

$$P_{\mathbf{\pi }}(\mathbf{\pi })={\mathrm{argmin}}_{{\mathbf{\pi }}^{\prime }}||\mathbf{\pi }-{\mathbf{\pi }}^{\prime }|{|}_2$$

$$\text{such that }{\mathbf{\pi }}^{\prime }⪰0,1^T{\mathbf{\pi }}^{\prime }=1.$$

$$\mathbf{\pi }\leftarrow P_{\mathbf{\pi }}\left(\mathbf{\pi }-\alpha \frac{\partial \ell }{\partial \mathbf{\pi }}\right),{\mathbf{\mu }}_c\leftarrow P_{{\mathbf{\mu }}_c}\left({\mathbf{\mu }}_c-\alpha \frac{\partial \ell }{\partial {\mathbf{\mu }}_c}\right),$$

$${\gamma }_1(\mathbf{x})=\frac{f(\mathbf{x}|{\mathbf{\mu }}_1)}{f(\mathbf{x}|{\mathbf{\mu }}_2)},{\gamma }_2(\mathbf{x})=1.$$

$$p^{\ast }(\mathbf{x})={\pi }_1^{\ast }\mathcal{N}(\mathbf{x}|{\mathbf{\mu }}^{\ast },I)+{\pi }_2^{\ast }\mathcal{N}(\mathbf{x}|-{\mathbf{\mu }}^{\ast },I)$$

$${\tilde{q}}_1(\mathbf{x})$$

$$Z_1={\pi }_1^{\ast }{e}^{{\mathbf{b}}^T({\mathbf{\mu }}^{\ast }-{\mathbf{\mu }}_2)}+{\pi }_2^{\ast }{e}^{-{\mathbf{b}}^T({\mathbf{\mu }}^{\ast }+{\mathbf{\mu }}_2)}.$$

$$Z_1={\pi }_1^{\ast }{e}^{2{\pi }_2^{\ast }{\mathbf{b}}^T{\mathbf{\mu }}^{\ast }}+{\pi }_2^{\ast }{e}^{-2{\pi }_1^{\ast }{\mathbf{b}}^T{\mathbf{\mu }}^{\ast }}.$$

$${\pi }_1^{\prime }={\pi }_1^{\ast }{Z}_1^{-1}{e}^{2{\pi }_2^{\ast }{\mathbf{b}}^T{\mathbf{\mu }}^{\ast }},{\pi }_2^{\prime }={\pi }_2^{\ast }{Z}_1^{-1}{e}^{-2{\pi }_1^{\ast }{\mathbf{b}}^T{\mathbf{\mu }}^{\ast }},$$

$$({\mathbf{b}}^{(t+1)}{)}^T{\mathbf{\mu }}^{\ast }=({\mathbf{b}}^{(t)}{)}^T{\mathbf{\mu }}^{\ast }+\delta ||{\mathbf{\mu }}^{\ast }|{|}_2^2>({\mathbf{b}}^{(t)}{)}^T{\mathbf{\mu }}^{\ast },$$

$$\frac{⟨{\mathbf{\mu }}_1^{(t+1)},{\mathbf{\mu }}^{\ast }⟩}{||{\mathbf{\mu }}_1^{(t+1)}||\cdot ||{\mathbf{\mu }}^{\ast }||}⩾\frac{⟨{\mathbf{\mu }}_1^{(t)},{\mathbf{\mu }}^{\ast }⟩}{||{\mathbf{\mu }}_1^t||\cdot ||{\mathbf{\mu }}^{\ast }||},$$

$$({\pi }_1+(\alpha /2)(Z_1-1),{\pi }_2-(\alpha /2)(Z_1-1)).$$

$${\mathbf{\mu }}_2\leftarrow {\mathbf{\mu }}_2+\alpha (\overline{\mathbf{x}}-{\mathbf{\mu }}_2).$$

$${\tilde{q}}_1(\mathbf{x})$$

$$Z_1$$

$$Z_1={\pi }_1^{\ast }\prod _{i=1}^D\left(1+{\pi }_2^{\ast }{\lambda }_i\right)+{\pi }_2^{\ast }\prod _{i=1}^D\left(1-{\pi }_1^{\ast }{\lambda }_i\right),$$

Full text

1 INTRODUCTION

The EM algorithm has a long history dating back to 1886 [1]. Its modern presentation was given in [2] and it has been applied widely on the maximum likelihood problem in mixture models, or equivalently, minimizing the negative log likelihood (also known as cross entropy). Some basic properties of EM have been derived. For example, EM monotonically decreases the loss [2] and it converges to critical points [3].

Mixture models [4] are important in clustering. The most popular one is the Gaussian mixture model (GMM). The optimization problem is non-convex, making the convergence analysis hard. This is true even in some simple settings such as one-dimensional or spherical covariance matrices. With separability assumptions, there has been a series of works proposing new algorithms for learning GMMs efficiently [5, 6, 7, 8, 9, 10, 11, 12, 13]. However, for the EM algorithm, not so much work has been done. One well-studied case is the two equally balanced mixtures [14, 15, 16] where local and global convergence guarantees are shown. For more than two components, [17] shows that there exist arbitrarily bad local minima. There has also been some study of EM on the convergence of Gaussian mixtures with more than two components [18, 19].

GMMs are well suited in practice for continuous variables. For discrete data clustering, we need discrete mixture models. One common assumption people make is naive Bayes, saying that in each component, the features are independent. This corresponds to the spherical covariance matrix setting in GMMs. Since discrete data can be converted to binary strings, in this paper we will focus on Bernoulli mixture models (BMMs), which have been applied in text classification [20] and digit recognition [21]. However, the theory is less studied. Although the Bernoulli and Gaussian distributions are both in the exponential family, we will show that BMMs have different properties from GMMs in the sense that GD is stable (i.e., gets trapped) at so-called one-cluster regions of BMMs, but unstable at (i.e., escapes) all one-cluster regions of GMMs. On the other hand, the two mixture models share similarities that we will explain.

In the practice of clustering with mixture models, it is common for algorithms to converge to regions with only $k$ clusters out of $m$ desired clusters (where $k<m$). We call them $k$-cluster regions. We observed in experiments that EM escapes such regions exponentially faster than GD for GMMs. For BMMs, GD may get stuck at some $k$-cluster regions and the probability can be very high, while EM always escapes such regions. These experimental results can be found in Section 6.

Theoretically, $k$-cluster regions are difficult to study. In this paper, we focus on one-cluster regions in mixtures of two components, which can be considered as a first step towards the more difficult problem. We succeed in explaining the different escape rates for GMMs with two components. We find that GD escapes one-cluster regions linearly while EM escapes exponentially fast. For BMMs, our theorem shows that there exist one-cluster regions where GD will converge to, given a small enough step size and a close enough neighborhood, but EM always escapes such regions exponentially. Such one-cluster regions can be arbitrarily worse than the global minimum. Experiments show that this contrast does not only hold for mixtures of two Bernoullis, but also BMMs with any number of components.

The rest of the paper is structured as follows. In Section 2, we give necessary background and notations. Following that, the main contributions are summarized in Section 3. In Section 4, we analyze the stability of EM and GD around one-cluster regions, for GMMs and BMMs with two components. The properties of such one-cluster regions are shown in Section 5. We provide supporting experiments in Section 6 and conclude in Section 7.

2 BACKGROUND AND NOTATIONS

A mixture model of $m$ components has the distribution:

[TABLE]

with the mixing vector ${\bm{\pi}}:=(\pi_{1},\dots,\pi_{m})$ on the $m-1$ dimensional probability simplex [22]: $\Delta_{m}:=\{\bm{\pi}:\bm{\pi}\succeq 0,\,{\bf 1}^{T}\bm{\pi}=1\}$. Here $f(\cdot|\boldsymbol{\mu}_{c})$ is the conditional distribution given cluster $c$, and parameters $\boldsymbol{\mu}_{c}$. The sample space (i.e., space of x) is $D$-dimensional.

We study the population likelihood, where there are infinitely many samples and the sample distribution is the true distribution. This assumption is common (e.g., in [14]), and it allows us to separate estimation error from optimization error. Since we consider the problem of minimizing the negative log likelihood, the loss function is:

[TABLE]

where the expectation is over $p^{*}(\mathbf{x})$, the true distribution. We assume by default this notation of expectation in this section. The loss function above is also known as cross entropy loss.

The true distribution $p^{*}(\mathbf{x})$ is assumed to have the same form as the model distribution: $p^{*}(\mathbf{x})=\sum_{c=1}^{m}\pi^{*}_{c}f({\mathbf{x}|\boldsymbol{\boldsymbol{\mu}}^{*}_{c}})$, where $\pi_{c}^{*}\in(0,1)$. Notice that the population assumption is not necessary in this section. We could replace $p^{*}(\mathbf{x})$ with a finite sample distribution, with i.i.d. samples drawn from $p^{*}(\mathbf{x})$.

Gaussian mixtures

For Gaussian mixtures, we consider the conditional distributions:

[TABLE]

where $\mathbf{x}\in\mathbb{R}^{D}$, $\boldsymbol{\mu}_{c}\in\mathbb{R}^{D}$ and the covariance matrix of each cluster ${\bm{\Sigma}}_{c}=I$.

Bernoulli mixtures

For Bernoulli mixtures, we have $\mathbf{x}\in\{0,1\}^{D}$, $\boldsymbol{\mu}_{c}\in[0,1]^{D}$ and:

[TABLE]

where the Bernoulli distribution is denoted as $B(x|\mu)=\mu^{x}(1-\mu)^{1-x}$. Slightly abusing the notation, we also use $B({\mathbf{x}|\boldsymbol{\mu}_{c}})$ to represent the joint distribution of $x_{1},\dots,x_{D}$ as a product of the marginal distributions. The expectation of the conditional distribution $B(\mathbf{x}|\boldsymbol{\mu}_{c})$ is:

[TABLE]

so $\boldsymbol{\mu}_{c}$ can be interpreted as the mean of cluster $c$. The covariance matrix of cluster $c$ is ${\bm{\Sigma}}_{c}={\rm diag}(\boldsymbol{\mu}_{c}*({\bf 1}-\boldsymbol{\mu}_{c}))$, where we use $*$ to denote element-wise multiplication. Also, we always assume $\boldsymbol{\mu}_{1}^{*}\in(0,1)^{D}$, $\boldsymbol{\mu}_{2}^{*}\in(0,1)^{D}$ and thus $\overline{{\bf x}}:=\mathds{E}[\mathbf{x}]\in(0,1)^{D}$.

2.1 EM algorithm

We briefly review the EM algorithm [4] for mixture models. Define the responsibility function111This is slightly different from the usual definition by a factor of $\pi_{c}$ (see, e.g., [4]). $\gamma_{c}(\mathbf{x})$:

[TABLE]

The EM update map $M$ is:

[TABLE]

Sometimes, it is better to interpret (2.7) in a different way. Define an unnormalized distribution $\tilde{q}_{c}(\mathbf{x})=p^{*}(\mathbf{x})\gamma_{c}(\mathbf{x})$. The corresponding partition function $Z_{c}$ and the normalized distribution $q_{c}(\mathbf{x})$ can be written as:

[TABLE]

where the integration is replaced with summation, given discrete mixture models. (2.7) can be rewritten as:

[TABLE]

This interpretation will be useful for our analysis near one-cluster regions in Section 4.

2.2 Gradient Descent

GD and its variants are the default algorithms in optimization of deep learning [23]. In mixture models, we have to deal with constrained optimization, and thus we consider projected gradient descent (PGD) [24]. From (2.1), (2.2) and (2.6), the derivative over $\bm{\pi}$ is:

[TABLE]

With (2.3) and (2.4), the derivatives over $\boldsymbol{\mu}_{c}$ are:

[TABLE]

The equation above is valid for both GMMs and BMMs. However, for BMMs, ${\bm{\Sigma}}_{c}^{-1}$ raises numerical instability near the boundary of $[0,1]^{D}$. Therefore, we use the following equivalent formula instead in the implementation:

[TABLE]

After a GD step, we have to project $(\bm{\pi},\boldsymbol{\mu}_{c})$ back into the feasible space $\Delta_{m}\times[0,1]^{D}$. For BMMs, we use Euclidean norm projection. The projection of the $\boldsymbol{\mu}_{c}$ part is simply:

[TABLE]

with $\min$ and $\max$ applied element-wise. Projection of $\boldsymbol{\mu}_{c}$ is not needed for GMMs. For the $\bm{\pi}$ part, we borrow the algorithm from [25], which essentially solves the following optimization problem:

[TABLE]

The projected gradient descent is therefore:

[TABLE]

with $\alpha$ the step size.

2.3 $k$-cluster region

We define a $k$-cluster region and a $k$-cluster point as:

Definition 2.1**.**

A $k$-cluster region is a subset of the parameter space where $||{\bm{\pi}}||_{0}=k$, with $||\cdot||_{0}$ denoting the number of non-zero elements. An element in a $k$-cluster region is called a $k$-cluster point.

In this work, we focus on one-cluster regions and mixtures of two components. In Section 6, we will show experiments on $k$-cluster regions for an arbitrary number of components. A key observation throughout our theoretical analysis and experiments, is that with random initialization, GD often converges to a $k$-cluster region where $k<m$, whereas EM almost always escapes such $k$-cluster regions, with random initialization and even with initialization in the neighborhood of $k$-cluster regions.

Now, let us study one-cluster regions for mixture models of two components. WLOG, assume $\pi_{1}=0$. To study the stability near such regions, we consider $\pi_{1}=\epsilon$, with $\epsilon$ sufficiently small. From (2.6), the responsibility functions can be approximated as:

[TABLE]

Under this approximation, $q_{2}(\mathbf{x})=p^{*}(\mathbf{x})$. For EM , $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ within one step based on (2.7). For GD, $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ at a linear rate222In the context of optimization, “a sequence $\{x_{k}\}_{k=1}^{\infty}$ converges to $L$ at a linear rate” means that there exists a number $\mu\in(0,1)$ such that $\lim_{k\to\infty}||x_{k+1}-L||/||x_{k}-L||=\mu$. based on (2.12)333For BMMs, assume $\boldsymbol{\mu}_{2}$ is not initialized on the boundary. The rate is upper bounded by the initialization as $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}\in(0,1)^{D}$..

For convenience, we define $\mathbf{b}:=\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}$ the difference between $\boldsymbol{\mu}_{1}$ and $\boldsymbol{\mu}_{2}$.

2.4 Stable fixed point and stable fixed region

A fixed point $p$ is stable under map $M$ if there exists a small enough feasible neighborhood $B$ of $p$ such that for any $p^{\prime}\in B$, $\lim_{n\to\infty}M^{n}(p^{\prime})=p$, where we use $M^{n}$ to denote a composition for $n$ times. Otherwise, the fixed point is called unstable. Similarly, a fixed region $R$ is stable under map $M$ if there exists a small enough feasible neighborhood $B$ of $R$ such that for any $p^{\prime}\in B$, $\lim_{n\to\infty}M^{n}(p^{\prime})$ exists and $\lim_{n\to\infty}M^{n}(p^{\prime})\in R$. We will use these definitions in Section 4.

3 MAIN CONTRIBUTIONS

We summarize the main contributions in this paper, starting from mixtures of two Gaussians. This is a well-known model that has been widely studied. However, we are not aware of any result regarding $k$-cluster points. As a first result in this line of research, we show that EM is better than GD in escaping one-cluster regions:

Result 3.1**.**

Consider a mixture of two Gaussians with $\pi_{1}^{*}\in(0,1)$, unit covariance matrices and true distribution

[TABLE]

When $\pi_{1}=\epsilon$ is initialized to be small enough, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$, EM increases $\pi_{1}$ exponentially fast, while GD increases $\pi_{1}$ linearly444Here, we use the usual notion of function growth. Do not confuse it with the rate of convergence..

This result indicates that EM escapes the neighborhood of one cluster-regions faster than GD by increasing the probability $\pi_{1}$ of cluster 1 at a rate that is exponentially faster than GD. The escape rates of EM and GD are formally proven in Theorems 4.1 and 4.2.

Our second result concerns mixtures of two Bernoullis:

Result 3.2**.**

For mixtures of two Bernoullis, there exist one-cluster regions that can trap GD, but do not trap EM. If any one-cluster region traps EM, it will also trap GD.

This result is formally stated in Theorems 4.3 and 4.7. It shows that EM is also better than GD in escaping one-cluster regions in BMMs. Our third result concerns the value of one-cluster regions, which is an informal summary of Theorems 5.1 and 5.2:

Result 3.3**.**

The one-cluster regions stated in Result 3.2 are local minima, and they can be $O(D)$ worse than the optimal value.

So far, we have seen that EM is better than GD in mixtures of two components, in the sense that EM escapes one-cluster regions exponentially faster than GD in GMMs, and that EM escapes local minima that trap GD in BMMs. Empirically, we show that this is true in general. In Section 6, we find that for BMMs with an arbitrary number of components and features, when we initialize the parameters randomly, EM always converges to an $m$-cluster point, i.e., all clusters are used in the model. Comparably, GD converges to a $k$-cluster point with a high probability, where only some of the clusters are employed in fitting the data (i.e., $k<m$). This result, combined with our analysis for mixtures of two components, implies that EM can be more robust than GD in terms of avoiding certain types of bad local optima when learning mixture models.

4 ANALYSIS NEAR ONE-CLUSTER REGIONS

In this section, we analyze the stability of EM and GD near one-cluster regions. Our theoretical results for EM and GD are verified empirically in Section 6.

4.1 Mixture of two Gaussians

We first consider mixtures of two Gaussians with identity covariance matrices, under the infinite-sample assumption. Due to translation invariance, we can choose the origin to be the midpoint of the two cluster means, such that: $\boldsymbol{\mu}_{1}^{*}=\boldsymbol{\mu}^{*}$ and $\boldsymbol{\mu}_{2}^{*}=-\boldsymbol{\mu}^{*}$.

4.1.1 EM algortihm

As long as $\pi_{1}$ is sufficiently small, the following theorem shows that EM will increase $\pi_{1}$ and therefore will not converge to a one-cluster solution. Recall from (2.9) that $\pi_{1}$ is multiplied by $Z_{1}$ at every step of EM and therefore we show that $Z_{1}>1$ almost everywhere, ensuring that $\pi_{1}$ will increase. Since $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ within one step, we only consider those points where $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$.

We can show that under mild assumptions, EM escapes one-cluster regions exponentially fast:

Theorem 4.1**.**

Consider a mixture of two Gaussians with $\pi_{1}^{*}\in(0,1)$, unit covariance matrices and true distribution

[TABLE]

When $\pi_{1}=\epsilon$ is initialized to be small enough, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$, EM increases $\pi_{1}$ exponentially fast.

Proof.

It suffices to show that $Z_{1}>1$ and $Z_{1}$ grows if initially $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$. To calculate $Z_{1}$, we first compute $\tilde{q}_{1}(\mathbf{x})=p^{*}(\mathbf{x})\gamma_{1}(\mathbf{x})$:

[TABLE]

This equation shows that $\tilde{q}_{1}(\mathbf{x})$ corresponds to an unnormalized mixture of Gaussians with their means shifted by $b$ and their mixing coefficients rescaled in comparison to $p^{*}(\mathbf{x})$. If $\mathbf{b}={\bf 0}$, then $\tilde{q}_{1}(\mathbf{x})=p^{*}(\mathbf{x})$. So, $\mathbf{b}$ describes how different $\tilde{q}_{1}(\mathbf{x})$ deviates from $p^{*}(\mathbf{x})$. The partition function can be computed by integrating out $\tilde{q}_{1}(\mathbf{x})$:

[TABLE]

In fact, $\boldsymbol{\mu}_{2}=\overline{\mathbf{x}}=(\pi_{1}^{*}-\pi_{2}^{*})\boldsymbol{\mu}^{*}$ which can be derived from $p^{*}(\mathbf{x})$. So, $Z_{1}$ becomes:

[TABLE]

Using the fact $e^{x}\geqslant 1+x$ and that equality holds iff $x=0$, we can show that when $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$, $Z_{1}>1$.

Now, let us show that $Z_{1}$ increases. It suffices to prove that $|\mathbf{b}^{T}\boldsymbol{\mu}^{*}|$ increases. From (4.1.1) and (2.9), we have the update equation for $\boldsymbol{\mu}_{1}$: $\boldsymbol{\mu}_{1}\leftarrow(\pi^{\prime}_{1}-\pi^{\prime}_{2})\boldsymbol{\mu}^{*}+\mathbf{b}$, with

[TABLE]

If $\mathbf{b}^{T}\boldsymbol{\mu}^{*}>0$, then $\pi^{\prime}_{1}>\pi_{1}^{*}$ and $\pi^{\prime}_{2}<\pi_{2}^{*}$. So, $\boldsymbol{\mu}_{1}^{(t+1)}=\boldsymbol{\mu}_{1}^{(t)}+\delta\boldsymbol{\mu}^{*}$ with $\delta=\pi^{\prime}_{1}-\pi_{1}^{*}+\pi_{2}^{*}-\pi^{\prime}_{2}>0$, and

[TABLE]

where we use the superscript $(t+1)$ to denote the updated values and $(t)$ to denote the old values.

Similarly, we can prove that $\mathbf{b}^{T}\boldsymbol{\mu}^{*}$ will decrease if $\mathbf{b}^{T}\boldsymbol{\mu}^{*}<0$. Hence, from (4.4), $Z_{1}$ increases under EM if initially $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$. ∎

This theorem above can be extended to any mixture of two Gaussians with known fixed covariance matrices $\Sigma$ for both clusters, using the transformation $\boldsymbol{\mu}_{i}\to\Sigma^{-1/2}\boldsymbol{\mu}_{i}$ and $\boldsymbol{\mu}^{*}\to\Sigma^{-1/2}\boldsymbol{\mu}^{*}$. We will prove it formally in Appendix C.

One may wonder what happens if originally $\mathbf{b}^{T}\boldsymbol{\mu}^{*}=0$. If we choose $\boldsymbol{\mu}_{1}$ randomly, $\mathbf{b}^{T}\boldsymbol{\mu}^{*}=0$ happens with probability zero since the corresponding Lebesgue measure is zero. Moreover, in numerical calculation such points are extremely unstable and thus unlikely. A similar condition appears in balanced mixtures of Gaussians as well, e.g., Theorem 1 in [15].

Rotation

The proof of Theorem 4.1 also shows an interesting phenomenon: $\boldsymbol{\mu}_{1}$ rotates towards $\boldsymbol{\mu}^{*}$ or $-\boldsymbol{\mu}^{*}$ (see Fig. 1). In fact, one can show that if $(\mathbf{b}^{(t)})^{T}\boldsymbol{\mu}^{*}>0$, then:

[TABLE]

where the equality holds iff $\boldsymbol{\mu}_{1}$ is a multiple of $\boldsymbol{\mu}^{*}$. A similar result can be shown for $(\mathbf{b}^{(t)})^{T}\boldsymbol{\mu}^{*}<0$.

Effect of Separation

Another interesting observation is that one can understand how the degree of separation affects the escape rate. If $||\boldsymbol{\mu}^{*}||$ is small, from (4.4), $Z_{1}$ will be close to $1$. Also, from (4.5) and (4.6), $\mathbf{b}^{T}\boldsymbol{\mu}^{*}$ will change more slowly, adding to the effect of slowing down the escape rate.

4.1.2 GD algorithm

Now, let us analyze the behavior of GD near one-cluster regions. After one GD iteration, the mixing coefficients become $(\pi_{1}+\alpha Z_{1},\pi_{2}+\alpha)$ before projection (where $\pi_{2}+\alpha Z_{2}=\pi_{2}+\alpha$ since $Z_{2}=1$). Assume a small step size $\alpha$ such that $\pi_{2}+\alpha>\pi_{1}+\alpha Z_{1}$. After projection based on (2.2), the updated mixing coefficients are:

[TABLE]

Combining with (2.9), whenever $Z_{1}>1$, both EM and GD increase $\pi_{1}$. In this sense, EM and GD achieve an agreement. However, in GD $\boldsymbol{\mu}_{1}$ changes little since the change is proportional to $\pi_{1}$. This argument is true for any mixture models with two components.

For GMMs, the update of $\boldsymbol{\mu}_{2}$ is:

[TABLE]

Hence, $\boldsymbol{\mu}_{2}^{(t)}$ is on the line segment between $\overline{{\bf x}}$ and $\boldsymbol{\mu}_{2}^{(0)}$, and $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ at a linear rate $1-\alpha$.

If $\boldsymbol{\mu}_{2}\neq\overline{{\bf x}}$, it is possible to have $Z_{1}<1$. For example, take $D=1$, $\mu_{2}=3$, $\mu_{1}=5$, $\pi_{1}^{*}=\pi_{2}^{*}=1/2$ and $\mu^{*}=2$, we have $Z_{1}=(1/2)(e^{-10}+e^{-2})<1$. In such cases, $\pi_{1}$ could either increase or decrease.

In the worst case, $\pi_{1}$ stays small until $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$. Then, from (4.4), $Z_{1}>1$ given $\mathbf{b}^{T}\boldsymbol{\mu}^{*}\neq 0$. Assuming $\boldsymbol{\mu}_{1}$ changes little in this process, we can conclude that with probability one GD escapes one-cluster regions. However, the growth of $\pi_{1}$ is linear compared to EM. Therefore, we have the following theorem:

Theorem 4.2**.**

For mixtures of two Gaussians and arbitrary initialization $(\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$, there exists a $\pi_{1}$ small enough such that when the mixture is initialized at $(\pi_{1},\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$, GD increases $\pi_{1}$ as a linear function of the number of time steps.

Proof.

For $\pi_{1}$ small enough, $\boldsymbol{\mu}_{1}$ stays in a small neighborhood, and $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ according to (4.9). Around $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$, $Z_{1}$ stays close to a constant and $Z_{1}\geqslant 1$. The linear growth of $\pi_{1}$ follows from (4.8). ∎

Notice that (4.8) is true for both GMMs and BMMs. The only difference is the computation of $Z_{1}$. From this equation, we know that if any one-cluster region traps EM, then $Z_{1}<1$ and $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$. On this configuration, for the GD algorithm, we know that $\pi_{1}$ decreases if it is not zero, and $\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}$ stays the same. Hence, this region traps GD as well. So, the following theorem holds:

Theorem 4.3**.**

For two-component mixture models, if any one-cluster region is stable for EM, it is stable for GD.

Proof.

Near any one-cluster region that is stable for EM, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and from (2.9), $Z_{1}<1$. Also, $\boldsymbol{\mu}_{1}$ changes very little since the update is proportional to $\pi_{1}$. So, $Z_{1}$ is a constant and $\pi_{1}$ decreases linearly according to (4.8). The parameters stop changing at $\pi_{1}=0$. ∎

4.2 Mixture of two Bernoullis

Now, let us see if Bernoulli mixtures have similar results. With the EM algorithm, $\boldsymbol{\mu}_{2}=\overline{\mathbf{x}}$ can be easily obtained. From the definition of $\tilde{q}_{c}(\mathbf{x})$, we have:

[TABLE]

Consider the unormalized binary distribution $B(x|\mu^{\prime}_{1},\mu^{\prime}_{2})=(\mu^{\prime}_{1})^{x}(\mu^{\prime}_{2})^{1-x}$ with $\mu^{\prime}_{1}+\mu^{\prime}_{2}\neq 1$. Then the normalization factor (partition function) is $\mu^{\prime}_{1}+\mu^{\prime}_{2}$. With this fact, one can derive that

[TABLE]

From the intuition of Gaussian mixtures, we similarly define $2{\bf\boldsymbol{\mu}}^{*}=\boldsymbol{\mu}_{1}^{*}-\boldsymbol{\mu}_{2}^{*}$ as the separation between the two true cluster means. As usual, ${\bf b}$ describes the separation between $\boldsymbol{\mu}_{1}$ and $\boldsymbol{\mu}_{2}$. With the notations above, we obtain:

[TABLE]

with

[TABLE]

Approximating $1+x\sim e^{x}$ and taking the covariance matrix to be the identity, we retrieve GMM (see (4.4)).

The update of ${\bm{\lambda}}$ can be computed from (2.9):

[TABLE]

where we denote $\Lambda_{i}=\mu_{1i}(1-\mu_{1i})$ and

[TABLE]

The derivation of (4.14) is presented in Appendix A.1. Notice that each $\lambda_{i}$ is an affine transformation of $\mu_{1i}$, as defined by (4.13). Hence, the domain of ${\bm{\lambda}}=(\lambda_{1},\dots,\lambda_{D})$ is a Cartesian product of closed intervals. We assume by default that ${\bm{\lambda}}$ is in the domain. Studying the update of ${\bm{\lambda}}$ instead of $\boldsymbol{\mu}_{1}$ is more convenient, as we will see in the following subsections that ${\bm{\lambda}}$ plays an important role in expressing the convergence results.

Also, we denote the covariance between features $i$ and $j$:

[TABLE]

We consider the case when there are no independent pairs, i.e., $\sigma_{ij}\neq 0$ for any $i\neq j$.

4.2.1 Attractive one-cluster regions for GD

When $Z_{1}<1$ and $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$, GD will get stuck in one-cluster regions, as shown in the proof of Theorem 4.3. If $\boldsymbol{\mu}_{2}\neq\overline{{\bf x}}$, $\boldsymbol{\mu}_{2}$ converges to $\overline{{\bf x}}$ at a linear rate:

[TABLE]

as can be seen from (2.12). Hence, we assume that $\boldsymbol{\mu}_{2}$ has already converged to $\overline{{\bf x}}$. With (4.12), we can define attractive one-cluster regions for GD:

Proposition 4.1** (attractive one-cluster regions for GD).**

Denote $T:=\{{\bm{\lambda}}|Z_{1}({\bm{\lambda}})<1\}$. The one-cluster regions $(\pi_{1},\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$ with ${\bm{\lambda}}(\boldsymbol{\mu}_{1})\in T$, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and $\pi_{1}=0$ are attractive for GD, given small enough step size. Here, we denote ${\bm{\lambda}}(\boldsymbol{\mu}_{1})$ as an affine function of $\boldsymbol{\mu}_{1}$, in the form of (4.13).

For example, when $D=2$, from (4.12), we have $Z_{1}=1+\pi_{1}^{*}\pi_{2}^{*}\lambda_{1}\lambda_{2}$. If $\lambda_{1}\lambda_{2}<0$, then $Z_{1}<1$, and GD will get stuck. If $\lambda_{1}\lambda_{2}>0$, $Z_{1}>1$, and then GD will escape one-cluster regions. In the latter case, ${\bm{\lambda}}\in\mathbb{R}_{++}^{2}$ or ${\bm{\lambda}}\in-\mathbb{R}_{++}^{2}$.

4.2.2 Positive regions

Now, let us study the behavior of EM. We will show that EM escapes one-cluster regions when ${\bm{\lambda}}\in\mathbb{R}_{++}^{2}\cup(-\mathbb{R}_{++}^{2})$, for the two-feature case. More generally with $D$ features, we should consider ${\bm{\lambda}}\in\mathbb{R}_{++}^{D}\cup(-\mathbb{R}_{++}^{D})$.

In the field of optimization, both $\mathbb{R}^{D}_{++}$ and $-\mathbb{R}^{D}_{++}$ are known as proper cones and a relevant order can be defined. They are critical in our analysis, so, we define these two cones as positive regions:

Definition 4.1** (positive regions).**

A positive region of ${\bm{\lambda}}$ is defined to be one of $\{\mathbb{R}^{D}_{++},-\mathbb{R}^{D}_{++}\}$. We denote it as $P:=\mathbb{R}^{D}_{++}\cup(-\mathbb{R}^{D}_{++})$.

The positive regions are interesting because $\mathbb{R}^{D}_{++}$ ensures that each element $\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}$ has the same sign as $\boldsymbol{\mu}_{1}^{*}-\boldsymbol{\mu}_{2}^{*}$ and $-\mathbb{R}^{D}_{++}$ ensures that each $\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{2}$ has the opposite sign of $\boldsymbol{\mu}_{1}^{*}-\boldsymbol{\mu}_{2}^{*}$.

Theorem 4.4**.**

For any number of features, $Z_{1}({\bm{\lambda}})>1$ in the positive regions.

Proof.

From (4.12),

[TABLE]

So, $\nabla Z_{1}({\bm{\lambda}})\succ 0$ if ${\bm{\lambda}}\succ 0$. Denote $g(t)=Z_{1}(t{\bm{\lambda}})$, then $Z_{1}({\bm{\lambda}})-1=g(1)-g(0)=\int_{0}^{1}{\bm{\lambda}}^{T}\nabla Z_{1}(t{\bm{\lambda}})dt>0$, given ${\bm{\lambda}}\succ 0$. we can similarly prove $Z_{1}({\bm{\lambda}})>1$ given ${\bm{\lambda}}\prec 0$. ∎

We also notice that positive regions are stable. Define $M({\bm{\lambda}})$ to be the EM update of ${\bm{\lambda}}$ based on the EM update of $\boldsymbol{\mu}_{1}$ and (4.13). For any number of features, the following lemma holds:

Lemma 4.1** (stability of positive regions).**

For any number of features, given $\pi_{1}=\epsilon$ small, the two positive regions $\mathbb{R}_{++}^{D}$ and $-\mathbb{R}_{++}^{D}$ are stable for EM, i.e., $M({\bm{\lambda}})\succ{\bm{\lambda}}$ for all ${\bm{\lambda}}\succ 0$ and $M({\bm{\lambda}})\prec{\bm{\lambda}}$ for all ${\bm{\lambda}}\prec 0$.

Proof.

If ${\bm{\lambda}}\succ 0$, then $B_{1i}>B_{2i}$ for any $i$ (based on (4.14) and (4.15)), and thus $M({\bm{\lambda}})\succ{\bm{\lambda}}$. A similar argument demonstrates that $M({\bm{\lambda}})\prec{\bm{\lambda}}$ when ${\bm{\lambda}}\prec 0$. ∎

From the lemma above, we conclude that in the positive regions, $|\mathbf{b}^{T}\boldsymbol{\mu}^{*}|$ will increase, similar to Theorem 4.1. Lemma 4.1 and Theorem 4.4 lead to the following:

Corollary 4.1**.**

For any number of features, if ${\bm{\lambda}}$ is initialized in a positive region, then EM will escape one-cluster regions exponentially fast.

Proof.

WLOG, we assume ${\bm{\lambda}}\succ 0$. From Theorem 4.4, at positive regions, $Z_{1}>1$. Define $g(t)=Z_{1}({\bm{\lambda}}+t(M({\bm{\lambda}})-{\bm{\lambda}}))$. $Z_{1}$ increases after an EM update because:

[TABLE]

where we use $M({\bm{\lambda}})\succ{\bm{\lambda}}$ from Lemma 4.1 and $\nabla Z_{1}({\bm{\lambda}})\succ 0$ if ${\bm{\lambda}}\succ 0$ from the proof of Theorem 4.4. From (2.9), $\pi_{1}$ increases exponentially. ∎

With GD, we have similar results as Corollary 4.1, but GD escapes one-cluster regions much more slowly, which can be derived in a similar way as in Section 4.1.

What if ${\bm{\lambda}}$ is not initialized in positive regions? The following theorem tells us that if $D=2$, no matter where ${\bm{\lambda}}$ is initialized, EM will almost always escape the one-cluster regions.

Theorem 4.5**.**

For $m=D=2$, given $\sigma_{12}\neq 0$ and $\overline{{\bf x}}\in(0,1)^{D}$, with the EM algorithm, $\pi_{1}=\epsilon$, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and uniform random initialization for $\boldsymbol{\mu}_{1}$, ${\bm{\lambda}}$ will converge to the positive regions at a linear rate with probability $1$. Therefore, EM will almost surely escape one-cluster regions.

Proof.

Here, we give a proof sketch, and the more detailed proof of Theorem 4.5 can be found in Appendix A.2. First, in the worst case, $||{\bm{\lambda}}||$ shrinks to a neighborhood of the origin. Then, near the origin, EM rotates ${\bm{\lambda}}$ towards positive regions. ∎

For general $m$ and $D$, we conjecture that EM almost always escapes $k$-cluster regions where $k<m$, as described in Appendix B.

4.2.3 EM as an ascent method

So far, we have studied positive regions $P$. We showed some nice properties of $P$ and proved that ${\bm{\lambda}}$ converges to $P$ almost everywhere when $D=2$. In this section, we show that EM can be treated as an ascent method for $Z_{1}$:

Theorem 4.6**.**

For ${\bm{\lambda}}$ in the feasible region, we have $\nabla Z_{1}({\bm{\lambda}})^{T}(M({\bm{\lambda}})-{\bm{\lambda}})\geqslant 0$, with equality holds iff ${\bm{\lambda}}=0$.

Proof.

From (4.12), $\nabla_{i}Z_{1}({\bm{\lambda}})=\pi_{1}^{*}\pi_{2}^{*}(B_{1i}-B_{2i})$. So, from (4.14), for any $i$, $\nabla_{i}Z_{1}({\bm{\lambda}})(M({\bm{\lambda}})_{i}-\lambda_{i})$ is nonnegative. Equality holds iff for all $i\in[D]$, $B_{1i}=B_{2i}$. It follows that for any $i$, $1+\pi_{1}^{*}\lambda_{i}=1-\pi_{2}^{*}\lambda_{i}$. Hence, equality holds iff ${\bm{\lambda}}=0$. ∎

With this theorem, we can take a small step in the EM update direction: ${\bm{\lambda}}+\alpha(M({\bm{\lambda}})-{\bm{\lambda}})$, with $\alpha$ small. In this way, we can always increase $Z_{1}$ until $Z_{1}>1$ and EM escapes the one-cluster regions.

4.3 EM vs GD

Now, we are ready to prove one of the main results in our paper: for mixtures of two Bernoullis, there exist one-cluster regions that can trap GD but not EM. We first prove a lemma similar to Lemma 4.1.

Lemma 4.2**.**

For all ${\bm{\lambda}}\neq 0$ and ${\bm{\lambda}}\succeq 0$, $M({\bm{\lambda}})\succ 0$ and $M({\bm{\lambda}})\succeq{\bm{\lambda}}$. Similarly, for all ${\bm{\lambda}}\neq 0$ and ${\bm{\lambda}}\preceq 0$, $M({\bm{\lambda}})\prec 0$ and $M({\bm{\lambda}})\preceq{\bm{\lambda}}$.

Proof.

WLOG, we prove the lemma for ${\bm{\lambda}}\succeq 0$. Notice that for any ${\bm{\lambda}}\succeq 0$, $M({\bm{\lambda}})\succeq{\bm{\lambda}}$ can be read directly from (4.14). For $\lambda_{i}=0$, $B_{1i}>B_{2i}$ since ${\bm{\lambda}}\neq 0$. Hence, $M({\bm{\lambda}})\succ 0$. ∎

This lemma leads to a main result in our work: for mixtures of two Bernoullis, there exist one cluster regions with nonzero measure that trap GD, but not EM.

Theorem 4.7**.**

For mixtures of two Bernoullis, given $\pi_{1}^{*}\in(0,1)$, $\overline{{\bf x}}\in(0,1)^{D}$ and $||\boldsymbol{\mu}^{*}||_{0}=D$, there exist one-cluster regions $B$ that trap GD, but EM escapes such one-cluster regions exponentially fast.

Proof.

Denote ${\bf e}_{i}$ as the $i^{\rm th}$ unit vector in the standard basis. For any $i\in[D]$, at ${\bm{\lambda}}={\lambda}_{i}{\bf e}_{i}$ and $\lambda_{i}>0$, from Lemma 4.2 and Theorem 4.4, we have $M({\bm{\lambda}})\succ 0$ and $Z_{1}(M({\bm{\lambda}}))>1$. Notice also that $Z_{1}({\bm{\lambda}})=1$ and $\nabla_{j}Z_{1}({\bm{\lambda}})>0$ for all $j\neq i$. Therefore, moving in the opposite direction of the gradient, it is possible to find a neighborhood $B$ of ${\bm{\lambda}}$, such that for all ${\bm{\lambda}}^{\prime}\in B$, $Z_{1}({\bm{\lambda}}^{\prime})<1$ and $Z_{1}(M({\bm{\lambda}}^{\prime}))>1$. This region traps GD due to Proposition 4.1, and EM escapes such one-cluster regions exponentially fast due to Corollary 4.1. ∎

5 ONE-CLUSTER LOCAL MINIMA

It is possible to show that some one-cluster regions are local minima, as shown in the following theorem:

Theorem 5.1**.**

The attractive one-cluster regions for GD defined in Proposition 4.1 are local minima of the cross entropy loss $\ell=-\mathds{E}[\log\left(\pi_{1}B(\mathbf{x}|\boldsymbol{\mu}_{1})+\pi_{2}B(\mathbf{x}|\boldsymbol{\mu}_{2})\right)]$.

Proof.

Impose $\pi_{2}=1-\pi_{1}$ and treat $\ell$ as a function of $(\pi_{1},\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$. Consider any small perturbation $(\delta\pi_{1},\delta\boldsymbol{\mu}_{1},\delta\boldsymbol{\mu}_{2})$. If $\delta\pi_{1}=0$, then the loss will not decrease as $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ is a local minimum of $-\mathds{E}[\log B(\mathbf{x}|\boldsymbol{\mu}_{2})]$. If $\delta\pi_{1}\neq 0$, the change of $\ell$ is determined by the first order. Only ${\partial\ell}/{\partial\pi_{1}}$ is nonzero, so the change is dominated by the first order:

[TABLE]

when $(\delta\pi_{1},\delta\boldsymbol{\mu}_{1},\delta\boldsymbol{\mu}_{2})$ is small enough. Therefore, we conclude that $(\pi_{1},\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$ is a local minimum. ∎

We call such minima one-cluster local minima. The following theorem shows that they are not global minima.

Theorem 5.2**.**

Assume $\boldsymbol{\mu}^{*}\neq 0$ and $\pi_{1}^{*}\in(0,1)$. For mixtures of two Bernoullis, one-cluster local minima cannot be global. The gap between the one-cluster local minima and the global minimum could be as large as ${\Theta}(D)$.

Proof.

The global minimum is obtained when $p(\mathbf{x})=p^{*}(\mathbf{x})$. Denote the optimal value as $\ell^{*}$. We have:

[TABLE]

with ${\rm KL}(p||q)$ the Kullback-Leibler divergence. Equality holds iff $p^{*}(\mathbf{x})=\prod_{i=1}^{D}p^{*}(x_{i})$, which means that the features are mutually independent. A direct consequence is that the features are pairwise independent. Nevertheless, we assumed at the end of Section 4.2 that $\sigma_{ij}\neq 0$ for any pair $(i,j)$. This is a contradiction.

Moreover, the difference could be very large. For example, taking $\pi_{1}^{*}=1/2$, $\boldsymbol{\mu}_{1}^{*}={\bf 1}$ and $\boldsymbol{\mu}_{2}^{*}=0$, the difference is $\ell^{*}-\ell_{1}=-(D-1)\log 2$. ∎

6 EXPERIMENTS

We analyzed the differences in behavior between EM and GD in Section 4 for mixtures of two components. For Gaussian mixtures with two components, we have done some experiments for the $D=2$ case with $10,000$ samples to approximate the population case, which can be shown in Figure 2. From this figure, we see that EM increases $\pi_{1}$ exponentially and GD increases $\pi_{1}$ linearly, when we initialize near a one-cluster region.

For Bernoulli mixtures, we aim to understand, at least empirically, how often EM and GD converge to a $k$-cluster point with $k<m$ with an arbitrary number of components. Since $k$-cluster points only utilize some of the components, this typically implies that the fixed points are bad. We will show that EM never converges to such bad critical points, whereas GD does.

We experiment with the infinite sample case for mixtures of Bernoullis of $m$ components with $D$ binary variables, where both $m$ and $D$ take values in $\{2,3,4,5,6\}$ for a total of 25 combinations. For each combination $(m,D)\in\{2,3,4,5,6\}^{2}$, We run (i) EM until convergence with at most $20000$ iterations; (ii) GD until convergence with a step size of $0.02$ and at most $10000$ iterations, for $60$ times each from random initialization. We find experimentally that EM always converges to an $m$-cluster point that is reasonably good compared to the optimal value (greater than $0.999$ in terms of the likelihood ratio). On the other hand, GD is much more likely to converge to a $k$-cluster point with $k<m$, especially as $m$ or $D$ increase. Table 1 summarizes our empirical findings for GD.

We also experiment with datasets of $2000$ samples generated from mixtures of Bernoullis of $m$ components with $D$ binary variables. We choose the same hyperparameters as in the previous experiment for EM and GD. As before, EM always converges to $m$-cluster points. The results for GD are summarized in Table 2. From this table, we see the similarity between infinite sample and finite sample cases. Hence, it is possible to extend our analysis for the population likelihood to a likelihood with a finite number of samples.

A remaining question is about the quality of the $k$-cluster points that GD converges to. Denote the loss of the $k$-cluster point that GD converges to as $\ell_{\rm GD}$, and the optimal loss as $\ell^{*}$. In Table 3, we compute $\exp(\ell^{*}-\ell_{\rm GD})$ as the ratio of the likelihood of the $k$-cluster points that GD converges to, to the optimal likelihood. We randomly choose 10 true distributions, compute the ratio and take the average, as shown in Table 3. We also take the worst-case ratio of the 10 true distributions, summarized in Table 4. From these two tables, we see that the $k$-cluster points become worse as $D$ increases. This agrees with our intuition for one-cluster points in Theorem 5.2. Also, with more clusters, the ratio is larger since it is easier to fit the data, especially when $D$ is small, and thus a $k$-cluster point may still behave well in this case.

7 CONCLUSIONS

In this paper, we identified $k$-cluster regions and studied one-cluster regions in mixture models of two components carefully. For mixtures of two Gaussians, when initialized around one-cluster regions, EM escapes such regions exponentially faster than GD, as also shown in experiments. For mixtures of two Bernoullis, we proved that there exist one-cluster local minima that always trap GD, but EM can escape such minima. For Bernoulli mixtures with a general number of components and any number of features, we showed experimentally that with random initialization, EM always converges to an $m$-cluster point where all components are used, but GD often converges to a $k$-cluster point where $k<m$. This means that only a subset of the components are employed.

Our work opens up a new direction of research. It would be interesting to know theoretically if EM almost always escapes $k$-cluster regions in mixture models. A conjecture for BMMs is given in Appendix B.

Acknowledgements

Guojun would like to thank Stephen Vavasis and Yaoliang Yu for useful discussions.

Appendix A Proofs for mixtures of Bernoullis

A.1 Derivation of (4.14)

From (4.2), the update of $\boldsymbol{\mu}_{1}$ is:

[TABLE]

where $F_{i}=\pi_{1}^{*}\mu_{1i}^{*}B_{1i}+\pi_{2}^{*}\mu_{2i}^{*}B_{2i}$ and $B_{1i},B_{2i}$ defined in (4.15). So,

[TABLE]

Bringing in the definition of $Z_{1}$ in (4.12), we have

[TABLE]

With the definitions of $\overline{x}_{i}$, $\mu_{i}^{*}$ and $\lambda_{i}$, we obtain:

[TABLE]

which, combined with (A.1), yields (4.14).

A.2 Proof of Theorem 4.5

In this section, we prove the following theorem:

Theorem 4.4**.**

For $m=D=2$, given $\sigma_{12}\neq 0$ and $\overline{{\bf x}}\in(0,1)^{D}$, with EM algorithm, $\pi_{1}=\epsilon$, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and uniform random initialization for $\boldsymbol{\mu}_{1}$, ${\bm{\lambda}}$ will converge to the positive regions at a linear rate with probability $1$. Therefore, EM will almost surely escape one-cluster regions.

We assume that $\sigma_{12}>0$ because the $\sigma_{12}<0$ can be similarly proved by relabeling $x_{2}\to 1-x_{2}$. The theorem is equivalent to showing that $\mathbf{b}$ converges to the regions where $b_{1}b_{2}>0$, due to (4.13) and Corollary 4.1. We also call these regions as positive regions. It is not hard to derive from (4.2) and (2.9) that the EM update is:

[TABLE]

with $\Lambda_{i}=\mu_{1i}(1-\mu_{1i})$.

We first notice some properties of $\sigma_{12}$, as can be easily seen from its definition. For convenience, in the following proof we define $\sigma:=\sigma_{12}S_{1}^{-1}S_{2}^{-1}$ which we call the normalized covariance.

Lemma A.1**.**

If $\sigma_{12}>0$, then $\sigma_{12}<\overline{x}_{1}(1-\overline{x}_{2})$, $\sigma_{12}<\overline{x}_{2}(1-\overline{x}_{1})$.

Proof.

Trivial from the definition of $\sigma_{12}$. ∎

A direct consequence is:

Corollary A.1**.**

If $\sigma_{12}>0$, then $\sigma^{2}S_{1}S_{2}<1$.

In the following lemma, we show that in a neighborhood of the origin, $\mathbf{b}$ almost always converges to the positive regions.

Lemma A.2** (Convergence with small $||\mathbf{b}||$ initialization, two features).**

Assume $\sigma_{12}>0$. $\exists\delta>0$ small enough, with a random $\mathbf{b}$ initialized from the $L_{1}$ ball $||\mathbf{b}||_{1}<\delta$ and the update function defined by EM, $\mathbf{b}$ converges to the positive regions $\{(b_{1},b_{2})|b_{1}b_{2}>0\}$ at a linear rate.

Proof.

In this case, $\Lambda_{i}$ and $Z$ are roughly constant, and

[TABLE]

The eigensystem of ${\bf A(\mathbf{b})}$ is:

[TABLE]

Applying ${\bf A}(\mathbf{b})$ for enough number of times, $\mathbf{b}$ will converge to a multiple of ${\bf v}_{1}$, where $\sigma b_{1}b_{2}>0$.

Let us make the argument above more concrete. Expand $\mathbf{b}$ as:

[TABLE]

Since $c_{1}=0$ is a measure zero set, we have $c_{1}\neq 0$ almost everywhere. WLOG, we assume $c_{1}>0$. Denote

[TABLE]

We first prove that each element of ${\bf A}(\mathbf{b})-{\bf A}(0)$ is bounded. This can be done by noticing:

[TABLE]

and that

[TABLE]

So, ${\bf A}(\mathbf{b})-{\bf A}(0)=O(\delta)$ in $||\mathbf{b}||_{1}<\delta$. From this fact, one can show in $||\mathbf{b}||_{1}<\delta$,

[TABLE]

with $c$ some constant. WLOG, we assume $\mathbf{b}$ is still in the negative region and thus $||\mathbf{b}||_{1}$ decreases, so our approximation is still valid. Applying EM for $k$ times and by use of (A.10), we know the result is at least (the generalized inequality is defined by the positive cone $\mathbb{R}^{2}_{++}$, see, e.g., [22])

[TABLE]

In the above analysis we used ${\bf A}(0){\bf u}\succeq 0$ for ${\bf u}\succeq 0$. For $\delta$ small enough, we know almost surely $c_{1}>{c\delta^{2}}/({\lambda_{1}-1})$. Therefore, at almost everywhere ${\bf A}(\mathbf{b})^{k}\mathbf{b}$ converges to the positive regions at a linear rate.

The choice of $L_{1}$ ball is irrelevant, since in finite vector space all $L_{p}$ norms are equivalent. ∎

Now, let us show that in the worst case $\mathbf{b}$ shrinks to a neighborhood of the origin. Hence, combined with Lemma A.2, we finish the proof. First, rewrite (A.4) and (A.5) as:

[TABLE]

The two contours $b^{\prime}_{1}=0,b^{\prime}_{2}=0$ are respectively:

[TABLE]

$f,g$ are both linear fractional functions of $b_{1}$, an example of which is depicted in Figure 3. The derivatives are:

[TABLE]

therefore, $f,g$ are both decreasing if $\sigma>0$. It follows that $f^{\prime}(0)=-(\sigma S_{1})^{-1}$ and $g^{\prime}(0)=-\sigma S_{2}$. From Corollary A.1,

[TABLE]

Now, let us look at the secant lines crossing the origin and $f(-\overline{x}_{1})$, $g(-\overline{x}_{1})$ separately. From (A.14) and (A.15),

[TABLE]

and one can obtain $f(1-\overline{x}_{1})$ and $g(1-\overline{x}_{1})$ similarly. So, $-b_{1}/b_{2}$ is bounded in the following region:

[TABLE]

and the bound is given by the slopes of the tangent lines at $(0,0)$ and the secant lines.

Another important point to notice is that $f(b_{1})$ and $g(b_{1})$ intersect exactly once. Which can be proved from Lemma A.1, (A.14) and (A.15):

Lemma A.3**.**

Assume $\sigma>0$, $f(b_{1})=g(b_{1})$ has exactly one solution $b_{1}=0$ in the feasible region $b_{1}\in[-\overline{x}_{1},1-\overline{x}_{1}]$.

Proof.

The solution $b_{1}=0$ is obvious. For $b_{1}\neq 0$, $f(b_{1})=g(b_{1})$ is equivalent to:

[TABLE]

We will show that the left hand side is always less than one. This is a linear function, so we only need to show it at both end points. At $b_{1}=-\overline{x}_{1}$, the left hand side can be simplified as:

[TABLE]

where we used Lemma A.1. Similarly, at $b_{1}=1-\overline{x}_{1}$, the left hand side of (A.21) is:

[TABLE]

∎

From Lemma A.3, the feasible region of $\mathbf{b}$ is divided into four parts by $f(b_{1})$ and $g(b_{1})$. If $f(b_{1})g(b_{1})>0$, then EM update goes to the positive region. Otherwise, $f(b_{1})g(b_{1})<0$. In this region, neither $b_{1}$ nor $b_{2}$ changes the sign. Because $-b_{1}/b_{2}$ is bounded, from (A.4) and (A.5), $||\mathbf{b}||^{(t+1)}_{1}\leqslant q||\mathbf{b}||^{(t)}_{1}$ with $0<q<1$ being a constant. So, in the worst case, $\mathbf{b}$ will converge to the neighborhood of the origin at a linear rate, and then shift to the positive regions at a linear rate, according to Lemma A.2. This lemma can be used because the EM update is not singular: it does not map a measure nonzero set to a measure zero set. Hence, after finitely many steps, at the neighborhood of the origin, the random distribution at the beginning is still random.

Figure 4 is an example of the trajectory.

Appendix B A General Conjecture

In this appendix, we propose a general conjecture for mixtures of Bernoullis:

Conjecture B.1**.**

For any number of clusters and a general number of features, with random initialization around $k$-cluster regions, EM will almost always converge to an $m$-cluster point.

Besides empirical evidence, we can show theoretical guarantees at $m=2$. In this case, the problem reduces to showing the convergence to positive regions proposed in Section 4.2.2. The convergence to the positive regions is observed empirically for mixtures of two Bernoullis, which always happens with random initialization.

The first result shows that if $||{\bm{\lambda}}||$ is small, ${\bm{\lambda}}$ will converge to the positive regions:

Proposition B.1** (Convergence with small $||{\bm{\lambda}}||$ initialization, general).**

For mixtures of two Bernoullis, $\exists\delta>0$ small enough, with a random ${\bm{\lambda}}$ initialized from the $L_{1}$ ball $||{\bm{\lambda}}||_{1}<\delta$ and the update function defined by EM, ${\bm{\lambda}}$ will converge to the positive regions.

Proof.

Around $||{\bm{\lambda}}||\sim 0$, $Z_{1}\sim 1$. Expanding (4.14) to the linear order, we find that the update of ${\bm{\lambda}}$ can be linearized as:

[TABLE]

where $A_{ii}=1$ and $A_{ij}=(2\mu_{i}^{*})^{2}\pi_{1}^{*}\pi_{2}^{*}S_{i}^{-1}Z_{1}^{-1}>0$ for $i\neq j$. Hence, ${\bf A}$ is an irreducible matrix. After enough iterations, ${\bm{\lambda}}$ will converge to the linear span of the largest eigenvector of ${\bf A}$.

From the Perron-Frobenius theorem [26], ${\bf A}$ has a unique largest real eigenvalue, and ${\rm eig}_{\rm max}({\bf A})\geqslant\min_{i}\sum_{j}{A_{ij}}>1$. Also, the maximal eigenvector of ${\bf A}$, ${\bf v}_{\rm max}$, is a multiple of an all positive vector. Therefore, we can prove the proposition in a similar fashion as the proof of Lemma A.2. ∎

This proposition also tells us that EM has an effect of rotating ${\bm{\lambda}}$ to the positive regions. It is interesting to observe that such unstable fixed point ${\bm{\lambda}}=0$ is analogous to the strict saddle points studied in [27]. It might be possible to use stable manifold theorem to prove our conjecture at $m=2$.

Another special case is when the $\max_{i}\lambda_{i}$ increases or $\min_{i}\lambda_{i}$ decreases. For a given ${\bm{\lambda}}$, we can order the components. WLOG, we assume $\lambda_{1}\leqslant\lambda_{2}\leqslant\dots\lambda_{i}<0=\lambda_{i+1}=\dots=\lambda_{j}<\lambda_{j+1}\leqslant\dots\leqslant\lambda_{D}$. We can show the following proposition:

Proposition B.2**.**

For mixtures of two Bernoullis, assume $\lambda_{1}\leqslant\lambda_{2}\leqslant\dots\lambda_{i}<0=\lambda_{i+1}=\dots=\lambda_{j}<\lambda_{j+1}\leqslant\dots\leqslant\lambda_{D}$, if $M({\bm{\lambda}})_{1}<\lambda_{1}$ or $M({\bm{\lambda}})_{D}>\lambda_{D}$, then ${\bm{\lambda}}$ will eventually converge to the positive regions. Otherwise, we have $M({\bm{\lambda}})_{1}\geqslant\lambda_{1}$ and $M({\bm{\lambda}})_{D}\leqslant\lambda_{D}$.

Proof.

From the definitions of $B_{1i}$ and $B_{2i}$, (4.15), we have

[TABLE]

$M_{1}({\bm{\lambda}})<\lambda_{1}$, from (4.14), tells us that $B_{11}-B_{21}<0$, and thus every $\lambda_{i}$ decreases. As each $\lambda_{i}$ decreases, $B_{1i}-B_{2i}$ will get smaller as well. Therefore, all $\lambda_{i}$’s decrease at least as a linear function. Since the feasible region is bounded, ${\bm{\lambda}}$ will converge to $-\mathbb{R}^{D}_{++}$ eventually.

Similarly, if $M_{D}({\bm{\lambda}})>\lambda_{D}$, we know that ${\bm{\lambda}}$ will converge to $\mathbb{R}^{D}_{++}$ eventually.

Otherwise, we must have $B_{11}-B_{21}\geqslant 0$ and $B_{1D}-B_{2D}\leqslant 0$, yielding $M_{1}({\bm{\lambda}})\geqslant\lambda_{1}$ and $M_{D}({\bm{\lambda}})\leqslant\lambda_{D}$. ∎

The two patterns $M({\bm{\lambda}})_{1}<\lambda_{1}$ and $M({\bm{\lambda}})_{D}>\lambda_{D}$ have been observed in experiments very frequently, while the case with $M({\bm{\lambda}})_{1}\geqslant\lambda_{1}$ and $M({\bm{\lambda}})_{D}\leqslant\lambda_{D}$ needs some further understanding.

Appendix C Proof for mixtures of two Gaussians with fixed equal covariances

In this appendix, we generalize our analysis for mixtures of two Gaussians with identity covariances to a slightly more general setting of mixtures of two Gaussians, given that the two covariance matrices are equal and fixed. The result is almost the same as Theorem 4.1, except that we replace the normal inner product in Euclidean space to a new inner product defined as $\langle a,b\rangle_{\Sigma}=a^{T}\Sigma^{-1}b$.

Theorem 4.4**.**

Consider a mixture of two Gaussians with $\pi_{1}^{*}\in(0,1)$, the same covariance matrices and true distribution

[TABLE]

When $\pi_{1}=\epsilon$ is initialized to be small enough, $\boldsymbol{\mu}_{2}=\overline{{\bf x}}$ and $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}\neq 0$, EM increases $\pi_{1}$ exponentially fast.

Proof.

It suffices to show that $Z_{1}>1$ and $Z_{1}$ grows if initially $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}\neq 0$. To calculate $Z_{1}$, we first compute $\tilde{q}_{1}(\mathbf{x})=p^{*}(\mathbf{x})\gamma_{1}(\mathbf{x})$:

[TABLE]

This equation shows that $\tilde{q}_{1}(\mathbf{x})$ corresponds to an un-normalized mixture of Gaussians with their means shifted by $b$ and their mixing coefficients rescaled in comparison to $p^{*}(\mathbf{x})$. If $\mathbf{b}={\bf 0}$, then $\tilde{q}_{1}(\mathbf{x})=p^{*}(\mathbf{x})$. So, $\mathbf{b}$ describes how different $\tilde{q}_{1}(\mathbf{x})$ deviates from $p^{*}(\mathbf{x})$.

The partition function can be computed by integrating out $\tilde{q}_{1}(\mathbf{x})$:

[TABLE]

In fact, $\boldsymbol{\mu}_{2}=\overline{\mathbf{x}}=(\pi_{1}^{*}-\pi_{2}^{*})\boldsymbol{\mu}^{*}$ which can be derived from $p^{*}(\mathbf{x})$. So, $Z_{1}$ becomes:

[TABLE]

Using the fact $e^{x}\geqslant 1+x$ and that equality holds iff $x=0$, we can show that when $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}\neq 0$, $Z_{1}>1$.

Now, let us show that $Z_{1}$ increases. It suffices to prove that $|\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}|$ increases. From (C), we have the update equation for $\boldsymbol{\mu}_{1}$: $\boldsymbol{\mu}_{1}\leftarrow(\pi^{\prime}_{1}-\pi^{\prime}_{2})\boldsymbol{\mu}^{*}+\mathbf{b}$, with

[TABLE]

If $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}>0$, then $\pi^{\prime}_{1}>\pi_{1}^{*}$ and $\pi^{\prime}_{2}<\pi_{2}^{*}$. So, $\boldsymbol{\mu}_{1}^{(t+1)}=\boldsymbol{\mu}_{1}^{(t)}+\delta\boldsymbol{\mu}^{*}$ with $\delta=\pi^{\prime}_{1}-\pi_{1}^{*}+\pi_{2}^{*}-\pi^{\prime}_{2}>0$, and

[TABLE]

where we use the superscript $(t+1)$ to denote the updated values and $(t)$ to denote the old values.

Similarly, we can prove that $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}$ will decrease if $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}<0$. Hence, from (C.4), $Z_{1}$ increases under EM if initially $\mathbf{b}^{T}\Sigma^{-1}\boldsymbol{\mu}^{*}\neq 0$. ∎