Maximum Number of Modes of Gaussian Mixtures

TL;DR

This paper investigates the maximum number of modes in Gaussian mixture models, providing new lower bounds and the first finite upper bound, advancing understanding of their complex density landscapes.

Contribution

It introduces improved lower bounds and the first finite upper bound on the number of modes in Gaussian mixtures, addressing a longstanding open problem.

Findings

Established new lower bounds on modes

Derived the first finite upper bound on modes

Enhanced understanding of Gaussian mixture density complexity

Abstract

Gaussian mixture models are widely used in Statistics. A fundamental aspect of these distributions is the study of the local maxima of the density, or modes. In particular, it is not known how many modes a mixture of $k$ Gaussians in $d$ dimensions can have. We give a brief account of this problem's history. Then, we give improved lower bounds and the first upper bound on the maximum number of modes, provided it is finite.

Equations

$$\varphi (x)=\frac{1}{\sqrt{\det (2\pi \Sigma )}}{\mathrm{e}}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}.$$

$${\Phi }_X(x)=\sum _{i=1}^k{\alpha }_i{\varphi }_{X_i}(x).$$

$$|{\mu }_2-{\mu }_1|\le 2\min ({\sigma }_1,{\sigma }_2)$$

$$|{\mu }_2-{\mu }_1|\le 2\sigma \sqrt{1+\frac{|\log (\alpha )-\log (1-\alpha )|}{2}}$$

$$x^{\ast }(\alpha )=[{\alpha }_1{\Sigma }_1^{-1}+{\alpha }_2{\Sigma }_2^{-1}+\dots +{\alpha }_k{\Sigma }_k^{-1}{]}^{-1}[{\alpha }_1{\Sigma }_1^{-1}{\mu }_1+{\alpha }_2{\Sigma }_2^{-1}{\mu }_2+\dots +{\alpha }_k{\Sigma }_k^{-1}{\mu }_k]$$

$$m(d,k)=\left(\genfrac{}{}{0ex}{}{d+k-1}{d}\right).$$

$${\Phi }_i(x):=\frac{1}{\sqrt{(2\pi {)}^d}\sqrt{{\delta }^3}}\exp \left(-\frac{1}{2{\delta }^3}|{\eta }_i(x){|}^2-\frac{1}{2}\Vert {\pi }_i(x)-{\mu }_i{\Vert }^2\right)$$

$$Q(p)=\left\{x\in {\mathbb{R}}^d\text{vline}|{\eta }_i(x)|\le \delta ,\forall 1\le i\le d\right\}$$

$$F_i^{\pm }=\left\{x\in {\mathbb{R}}^d\text{vline}{\eta }_i(x)=\pm \delta ,|{\eta }_j(x)|\le \delta \forall j\ne i\right\}.$$

$$\sqrt{{\delta }^3}\underset{x\in F_j^{\pm }}{\max }{\Phi }_i(x)$$

$$\sqrt{{\delta }^3}\underset{x\in F_i^{\pm }}{\max }{\Phi }_i(x)$$

$$\sqrt{{\delta }^3}\underset{x\in Q}{\max }{\Phi }_i(x)$$

$${\varphi }_i:=\frac{1}{\sqrt{(2\pi {)}^d}}\exp \left(-\frac{1}{2}\Vert p-{\mu }_i{\Vert }^2\right).$$

$$\sqrt{{\delta }^3}{\Phi }_i(x)\le \frac{1}{\sqrt{(2\pi {)}^d}}\exp \left(-\frac{1}{2}\Vert {\pi }_i(x)-{\mu }_i{\Vert }^2\right).$$

$$\sqrt{{\delta }^3}{\Phi }_i(x)\le \frac{1}{\sqrt{(2\pi {)}^d}}\exp \left(-\frac{1}{2{\delta }^3}{\delta }^2\right)\stackrel{\delta \to 0}{\to }0.$$

$$\sqrt{{\delta }^3}{\Phi }_i(x)\le \frac{1}{\sqrt{(2\pi {)}^d}}\exp \left(-\frac{1}{8{\delta }^3}|{\eta }_i(p){|}^2\right)\stackrel{\delta \to 0}{\to }0.$$

$$\sqrt{{\delta }^3}\underset{x\in \partial Q}{\max }\Phi (x)$$

$$\sqrt{{\delta }^3}\Phi (p)$$

$$2^{d+\left(\genfrac{}{}{0ex}{}{k}{2}\right)}(5+3d{)}^k.$$

$$d_1\cdots d_n(5+n+d_1+\dots +d_n{)}^k\cdot 2^{\frac{k(k-1)}{2}}.$$

$$\widehat{g_i}(x,t):=F_i(x_1,\dots ,x_n,y_1(x),\dots ,y_{k-1}(x),t{y}_k(x))$$

$$v_r:=(-1{)}^{n+1-r}\det \frac{\partial (\widehat{g_1},\dots ,\widehat{g_n})}{\partial (x_1,\dots ,x_{r-1},x_{r+1},\dots ,x_n,t)}$$

$$\widehat{g_1}=0,\dots ,\widehat{g_n}=0,v_{n+1}=\det \frac{\partial \hat{g}}{\partial x}=0.$$

$$N\le d_1\cdots d_n(n+D)(5+(n+1)+(D+n+D){)}^{k-1}\cdot 2^{\frac{(k-1)(k-2)}{2}}$$

$$\widehat{g_1}=0,\dots ,\widehat{g_n}=0,{\lambda }_1{x}_1+\dots +{\lambda }_n{x}_n+{\lambda }_{n+1}u+\mu =0,$$

$$q\le d_1\cdots d_n\cdot 1\cdot (5+(n+1)+(D+1){)}^{k-1}\cdot 2^{\frac{(k-1)(k-2)}{2}}$$

$$\begin{array}[]{rcl}N+q&\leq&d_{1}\cdots d_{n}\left[(n+D)(6+2n+2D)^{k-1}+(7+n+D)^{k-1}\right]\cdot 2^{\frac{(k-1)(k-2)}{2}}\\ &<&d_{1}\cdots d_{n}\left[(n+D)(5+n+D)^{k-1}\cdot 2^{k-1}+5(5+n+D)^{k-1}\cdot 2^{k-1}\right]\cdot 2^{\frac{(k-1)(k-2)}{2}}\\ &=&d_{1}\cdots d_{n}(5+n+D)^{k}\cdot 2^{\frac{k(k-1)}{2}}.\end{array}$$

$$g_0(x_0,x_1,\dots ,x_n)=x_0^2+x_1^2+\dots +x_n^2-R^2=0$$

$$N^{\prime }\le 2\cdot d_1\cdots d_n\cdot (n+1+D)(5+(n+2)+(2+D+n+1+D){)}^{k-1}\cdot 2^{\frac{(k-1)(k-2)}{2}}$$

$$q^{\prime }\le 2\cdot d_1\cdots d_n\cdot 1\cdot (5+(n+2)+2+D+1{)}^{k-1}\cdot 2^{\frac{(k-1)(k-2)}{2}}$$

$$\begin{array}[]{rcl}\dfrac{N^{\prime}+q^{\prime}}{2}&\leq&d_{1}\cdots d_{n}\left[(1+n+D)(10+2n+2D)^{k-1}+(10+n+D)^{k-1}\right]\cdot 2^{\frac{(k-1)(k-2)}{2}}\\ &<&d_{1}\cdots d_{n}\left[(1+n+D)(5+n+D)^{k-1}\cdot 2^{k-1}+4(5+n+D)^{k-1}\cdot 2^{k-1}\right]\cdot 2^{\frac{(k-1)(k-2)}{2}}\\ &=&d_{1}\cdots d_{n}(5+n+D)^{k}\cdot 2^{\frac{k(k-1)}{2}},\end{array}$$

$$2\cdots 2(5+d+2+\dots +2{)}^k\cdot 2^{\frac{k(k-1)}{2}}=2^d(5+d+2d{)}^k{2}^{\left(\genfrac{}{}{0ex}{}{k}{2}\right)}=2^{d+\left(\genfrac{}{}{0ex}{}{k}{2}\right)}(5+3d{)}^k.$$

Full text

1. Introduction

The $d$-dimensional Gaussian distribution $N(\mu,\Sigma)$ can be defined by its probability density function

[TABLE]

where $\mu\in\mathbb{R}^{d}$ is the mean vector and the symmetric positive definite $d\times d$ matrix $\Sigma$ is the covariance matrix. This density has a unique absolute maximum at $x=\mu$. Now consider a mixture distribution $X$ consisting of $k$ Gaussian components $X_{i}\sim N(\mu_{i},\Sigma_{i})$ and mixture weights $\alpha_{i}$ for $i=1,\ldots,k$, so that:

[TABLE]

This is again a probability density function given that $\alpha_{i}\geq 0$ and $\alpha_{1}+\alpha_{2}+\cdots+\alpha_{k}=1$. In other words, the density of a Gaussian mixture is a convex combination of Gaussian densities. Such mixtures can exhibit quite complex behavior even for a small number $k$ of components. This is a feature that makes them attractive for modeling in applications.

A fundamental property of a probability density function is the number of modes, i.e. local maxima, that it possesses. For Gaussian mixtures, this is especially relevant in applications such as clustering [9, p. 383]. For example, the mean shift algorithm converges if there are only finitely many critical points [17].

We will be interested in the maximal number $m(d,k)$ of local maxima for $d$-dimensional Gaussian mixtures with $k$ components. Shockingly, it is not known whether this maximal number is always finite for general Gaussian mixtures. On the other hand, we stress that the number of modes is a property of a Gaussian mixture density with fixed parameters and no sample involved; it should not be confused with the number of local maxima of the likelihood function of a Gaussian mixture model (a relevant but different question, see [2] and [10]).

Remark 1**.**

Since a single Gaussian has a unique global maximum at its mean $\mu$, we have that $m(d,1)=1$ for all $d\geq 1$.

Our goal is to look for lower and upper bounds for $m(d,k)$. Our main results in this direction are Theorem 6 for the former and Corollary 11 for the latter.

2. Background

As stated in the introduction, a single Gaussian has a unique mode, that is, $m(d,1)=1$ for all $d\geq 1$. The simplest case when there is an actual mixture has $d=1$ and $k=2$: a mixture of two univariate Gaussians $X_{1}\sim N(\mu_{1},\sigma_{1}^{2})$ and $X_{2}\sim N(\mu_{2},\sigma_{2}^{2})$, with mixture parameter $\alpha\in(0,1)$. It was observed historically that in this scenario the number of modes was either $1$ or $2$, with the following heuristics:

• •

If the distance between the component means is small, then the mixture is unimodal (independently of $\alpha$).

• •

If the distance between the component means is large enough, then there is bimodality unless $\alpha$ is close to [math] or $1$.

A.C. Cohen (1953) and Eisenberger (1964) obtained some first explicit conditions in these directions [3]. Notably, if $\alpha=\frac{1}{2}$ and $\sigma_{1}=\sigma_{2}=\sigma$, then the mixture is unimodal (with mode at $\frac{\mu_{1}+\mu_{2}}{2}$) if and only if $|\mu_{2}-\mu_{1}|\leq 2\sigma$. A few years later, J. Behboodian gave a proof that indeed $m(1,2)=2$ by showing the number of critical points of the density is at most three, and finds that

[TABLE]

is a sufficient condition for unimodality. Furthermore, if $\sigma_{1}=\sigma_{2}=\sigma$, then

[TABLE]

is again a sufficient condition for having only one mode in the mixture. Starting the 21st century, it was Carreira-Perpiñán and Williams who had particular interest in the problem [5]. Using scale-space theory, they prove that $m(1,k)=k$; any univariate Gaussian mixture with $k$ components has at most $k$ modes. A natural conjecture could be that $m(d,k)=k$ for all $d,k$, that is, a mixture with $k$ Gaussian components can have at most $k$ modes. 111In June 2016, a discussion thread on the ANZstat mailing list (e-mail bulletin board for statistics in Australia and New Zealand) with the title “an interesting counter-intuitive fact” referred to the fact that a Gaussian mixture can have more modes than components. However, this fails already when $d=k=2$, since a mixture of two bivariate Gaussians can have three distinct modes (and actually, $m(2,2)=3$).

Example 2**.**

Consider $X_{1}\sim N\left(\begin{pmatrix}1\\ 0\end{pmatrix},\begin{pmatrix}1&0\\ 0&0.1\end{pmatrix}\right)$ and $X_{2}\sim N\left(\begin{pmatrix}0\\ 1\end{pmatrix},\begin{pmatrix}0.1&0\\ 0&1\end{pmatrix}\right)$, with $\alpha=\frac{1}{2}$. There are two modes close to the original means at $(1,0)$ and $(0,1)$ and there is also a third mode near the origin. This situation is illustrated in the contour plot of Figure 2, with the two means marked ‘+’ and the three modes marked in red.

Special attention can be paid to assumptions on the variances. A mixture is said to be homoscedastic if all the variances in the components are equal: $\Sigma_{i}=\Sigma$ for $i=1,\ldots,k$. On the other hand, a mixture is said to be isotropic if $\Sigma_{i}=\sigma_{i}^{2}I$ for $i=1,\ldots,k$, so that covariances are scalar matrices and the densities have a ‘spherical’ shape. Note that, up to coordinate change, homoscedastic mixtures are (homoscedastic) isotropic. Carreira-Perpiñán and Williams conjectured in [5] that if one is restricted to homoscedastic Gaussian mixtures, then the maximum number of modes is actually $k$, and verified this numerically for many examples in a brute force search. Denoting by $h(d,k)$ this maximum number, they asserted that $h(d,k)=k$ for any $d,k\geq 1$.

Remark 3**.**

It holds that $h(d,k)\leq m(d,k)$ for all $d,k$ , and $h(d,k)$ is also the maximum number of possible modes of a Gaussian mixture with all unit covariances (by the note above and the fact that the number of modes remains invariant under affine transformations)

However, later J.J. Duistermaat emailed the authors of [6] with a counterexample in dimension $d=2$ with $k=3$ isotropic components, each on the vertex of an equilateral triangle. This configuration gives 4 modes for a small window of parameters, disproving the conjecture.

Example 4**.**

Consider the isotropic mixture with components $X_{1}\sim N((1,0),\sigma^{2}I_{2})$, $X_{2}\sim N((-\frac{1}{2},\frac{\sqrt{3}}{2}),\sigma^{2}I_{2})$ and $X_{3}\sim N((-\frac{1}{2},-\frac{\sqrt{3}}{2}),\sigma^{2}I_{2})$ with $\sigma^{2}=0.53$ and $\alpha=(\frac{1}{3},\frac{1}{3},\frac{1}{3})$. There is a mode on each of the three line segments between the origin and the means, and there is also a fourth mode at the origin. This situation is illustrated in the contour plot of Figure 3, with the three means marked ‘+’ and the four modes marked as red points.

In terms of contribution to the study of the topography of Gaussian mixture densities, S. Ray and B. Lindsay initiate a systematic study and ask interesting questions in [13]. They consider the ridgeline function $x^{*}:\Delta_{k}\rightarrow\mathbb{R}^{d}$ given by

[TABLE]

where $\Delta_{k}=\{(\alpha_{1},\alpha_{2},\dots,\alpha_{k})|\,\alpha_{i}\geq 0\text{ and }\alpha_{1}+\alpha_{2}+\dots+\alpha_{k}=1\}$ denotes the $(k-1)$-dimensional probability simplex, obtaining as its image the ridgeline variety $\mathcal{M}=\operatorname{Im}(x^{*})$ that contains all critical points of $\Phi_{X}$ for fixed $\mu_{1},\ldots,\mu_{k}$ and fixed $\Sigma_{1},\ldots,\Sigma_{k}$. This fact is useful, for example, in the case of homoscedastic mixtures, whose critical points (and in particular all modes) lie in the convex hull of the component means (a result that appeared first in [5]). It would be interesting to study the locus of critical points of the Gaussian mixture density function as the means and/or covariances vary.

In the conclusion of [13], the following line appears: “one might ask if there exists an upper bound for the number of modes, one that can be described as a function of $k$ and $d$”.

Assuming this bound is finite, we answer this question in the affirmative in Section 5.

3. Examples and Conjecture

The appearance of a possible extra mode in dimension $d=2$ when having $k=2$ components carries over to higher dimensions. Ray and Ren proved in [14] that $m(d,2)=d+1$. That is, one can get as many as $d+1$ modes from just a two component Gaussian mixture in dimension $d$. Looking for further progress, Ray proposed the maximum number of modes problem for the 2011 AIM Workshop on Singular Learning Theory, organized by Steele, Sturmfels and Watanabe [16]. The problem was discussed, and it led to the following conjecture:

Conjecture 5**.**

(Sturmfels, AIM 2011) For all $d,k\geq 1$,

[TABLE]

This conjecture matches correctly all the known values for $m(d,k)$ so far, which we have presented. In the next section we will show that for $d=2$ there exist Gaussian mixtures that achieve as many as $\binom{k+1}{2}$ modes, showing that (4) is a lower bound on $m(2,k)$.

Let $d=2$ and $k=3$. The conjectured bound gives $\binom{2+3-1}{2}=6$ modes. In Figure 4 we give an example of a Gaussian mixture that has this number of modes. The configuration relies on the deformation of 3 lines arranged in an equilateral triangle, and taking means as the middle points on the 3 sides, with all weights $\frac{1}{3}$. Apart from the modes coming from the means, the other 3 modes lie near the corresponding triangle vertices.

One could ask if in this anisotropic case, there exists a counterexample to Conjecture 5 in the spirit of Duistermaat’s mixture. Specifically, could an extra mode be formed at the origin for some values of the covariance parameters? This would give a total of 7 modes. Note that by rotational symmetry, the origin is always a critical point. Indeed, if the Gaussians are very concentrated on the lines (like in Figure 4), then the origin is a local minimum. If they diffuse enough, then it will eventually become a mode. The problem is that in this diffusion, the modes coming from the means quickly become saddle points, as illustrated in Figure 5. We argue that an intermediate scenario of a total of 7 modes is actually impossible. Indeed, consider any height of the equilateral triangle. Again by symmetry, the corresponding modes near the vertex and middle point of the triangle lie on this height. Restricting the Gaussian mixture density to that line, the components corresponding to the opposite sides project to the same kernel; thus obtaining a combination of two Gaussian kernels. Since we know that the number of modes is at most two in one dimension, not all three of the critical points lying on the line can be modes.

In [7], a construction of an isotropic (and homoscedastic) mixture of Gaussians is presented. One considers products of triangles, using Duistermaat’s counterexample with $4$ modes of Section 2 as the basic building block to obtain $4^{n}$ modes in dimension $d=2n$ with $k=3^{n}$ components. This gives an example where the number of modes is superlinear $k^{1.261}$ in the number of components (however, note that the dimension $d=2\log_{3}k$ also grows with $k$).

In the following section, we provide configurations for any choice of $d>1$ and $k>1$ having $\binom{k}{d}+k$ modes. If we let $d$ grow logarithmically with $k$ as in [7], we obtain superpolynomially (but subexponentially) many modes.

4. Many Modes

In this section, we prove that Gaussian mixtures can have many modes.

Theorem 6**.**

Given integers $k,d\geq 2$, there is a mixture of $k$ Gaussians in $\mathbb{R}^{d}$ with at least $\binom{k}{d}+k$ modes. That is, $m(d,k)\geq\binom{k}{d}+k$.

These are the terms $i=1$ and $i=d$ in the expansion $\binom{d+k-1}{d}=\sum_{i=1}^{d}\binom{d-1}{d-i}\binom{k}{i}$ of the conjectured bound (4). For $d=2$, our bound agrees with (4).

Proof.

Starting from a generic arrangement $H_{1},\ldots,H_{k}$ of $k$ affine hyperplanes in $\mathbb{R}^{d}$, we are going to define a family $\Phi=\Phi^{\delta}$ of Gaussian mixtures depending on a parameter $\delta>0$. Around each of the $\binom{k}{d}$ intersection vertices $p$ of the arrangement, we construct neighborhoods $Q=Q(p)$, also depending on $\delta$, so that for $\delta$ small enough, we have $\Phi^{\delta}|_{\partial Q}<\Phi^{\delta}(p)$. This certifies the existence of a mode in $Q$ for each $p$ (for this, we may assume that $k\geq d$ since otherwise there are no vertices). In addition, there will be a mode near each of the $k$ means. For each $i=1,\ldots,k$, denote by $\pi_{i}\colon\mathbb{R}^{d}\to H_{i}$ the orthogonal projection, and pick an affine map $\eta_{i}\colon\mathbb{R}^{d}\to\mathbb{R}$ such that $|\eta_{i}(x)|$ is the distance from $x$ to $H_{i}$. Further, choose means $\mu_{i}\in H_{i}$ outside the other $H_{j}$. Then, our $i$th component will be a standard Gaussian with mean $\mu_{i}$ along $H_{i}$ with variance $\delta^{3}$ in the direction normal to $H_{i}$:

[TABLE]

For the mixture, we take all coefficients to be equal: $\Phi=\frac{1}{k}\sum_{i}\Phi_{i}$. Let $p$ be one of the intersection vertices; without loss of generality, $\left\{p\right\}=H_{1}\cap\cdots\cap H_{d}$. For $\delta>0$, we define the neighborhood $Q$ of $p$ to be

[TABLE]

(note that $\eta_{i}(p)=0$ for $i=1,\dots,d$). This is an affine cube with center $p$. Now we consider each of its $2d$ facets $F_{i}^{\pm}$, $1\leq i\leq d$, where

[TABLE]

We will show that, around the point $p$, as $\delta\to 0$,

[TABLE]

where $\phi_{i}$ is the positive number

[TABLE]

To establish (6), observe that along $F_{j}^{\pm}$ ($j\neq i$), we have

[TABLE]

with equality at the center of $F_{j}^{\pm}$. The right hand side is a continuous function of $x$, it is independent of $\delta$, and it evaluates to $\phi_{i}$ at $p$. As all of $F_{j}^{\pm}$ converges to $p$, we must have (6).

To establish (7), observe that along $F_{i}^{\pm}$, we have

[TABLE]

To establish (8), fix $i\in\{d+1,\ldots,k\}$. Observe that $|\eta_{i}(p)|>0$. As the diameter of $Q$ is linear in $\delta$, for small enough $\delta$, we have $|\eta_{i}(x)|>\frac{1}{2}|\eta_{i}(p)|$ for all $x\in Q$. Hence, for those $\delta$’s

[TABLE]

Adding up (6), (7) and (8), we get for the mixture density $\Phi=\frac{1}{k}\sum_{i=1}^{k}\Phi_{i}$ that

[TABLE]

As the limit (9) is smaller (by one $\phi_{j}>0$) than the limit (10), there must be some $\delta^{\star}(p)>0$ so that for $0<\delta<\delta^{\star}(p)$ we have $\max_{x\in\partial Q}\Phi(x)<\Phi(p)$. Then the point $p^{\prime}$ where the continuous function $\Phi$ takes it maximum over the compact set $Q$ will be in the interior of $Q$, and hence is a local maximum. Choosing $\delta^{\star}$ to be the minimum over the $\delta^{\star}(p)$ over all intersection vertices $p$, we obtain a mixture with at least $\binom{k}{d}$ modes.

The argument for the existence of a mode near $\mu_{i}$ is similar, but much simpler. Fix a compact neighborhood $Q$ of $\mu_{i}$ which avoids the hyperplanes $H_{j}$ for $j\neq i$. For small $\delta$, as in (8), the $\Phi_{j}$ for $j\neq i$ become negligible along $Q$. Since the remaining $\Phi_{i}$ attains its maximum at $\mu_{i}$, the value of the mixture density $\Phi$ at $\mu_{i}$ will be larger than its values along $\partial Q$. Thus we obtain the existence of another $k$ modes for sufficiently small $\delta$. ∎

5. Not Too Many Modes

The main result of this section is to present an upper bound on the number of modes of a Gaussian mixture. We start by looking at the set of critical points and we will use Khovanskii’s theory on fewnomials, see [11].

Theorem 7**.**

For all $d,k\geq 1$, the number of non-degenerate critical points for the density of a mixture of $k$ Gaussians in $\mathbb{R}^{d}$ is bounded by

[TABLE]

This will follow from a Khovanskii-type theorem that bounds the number of nondegenerate solutions to a system of polynomial equations that includes transcendental functions. Such a version where the transcendental functions are exponentials of linear forms was first presented by Khovanskii to illustrate his theory of fewnomials [11, p.12]. In our case, however, we will be interested in exponentials of quadratic forms.

Theorem 8**.**

For $1\leq i\leq n$, let $F_{i}\in\mathbb{R}[x_{1},\ldots,x_{n},y_{1},\ldots,y_{k}]$ be polynomials of degree $d_{i}$ and for $1\leq j\leq k$ consider the exponential quadratic forms $y_{j}(x)=\mathrm{e}^{(x-\mu_{j})^{T}Q_{j}(x-\mu_{j})}$, with $\mu_{j}\in\mathbb{R}^{n}$ and $Q_{j}\in\mathbb{R}^{n\times n}$. If $g_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ are given by $g_{i}(x)=F_{i}(x_{1},\ldots,x_{n},y_{1}(x),\ldots,y_{k}(x))$ then the number of non-degenerate solutions to the system $g_{1}=g_{2}=\cdots=g_{n}=0$ is finite and bounded by

[TABLE]

In order to prove his theorem, Khovanskii gives first a sketch making simplifying assumptions (and skipping technical details) and fills the theory in his next two chapters. This sketch is also presented in [15] and [4], and we will present the proof of our theorem in the same way. We will need the following lemma (for a proof see e.g. Theorem 4.3 in [15]).

Lemma 9** (Khovanskii-Rolle).**

Let $C\subset\mathbb{R}^{n+1}$ be a smooth curve that intersects the hyperplane $H$ given by ${x_{n+1}=0}$ transversally, and $v=(v_{1},\ldots,v_{n+1}):C\rightarrow\mathbb{R}^{n+1}$ a smooth nonvanishing tangential vector field to $C$. Then $|(C\cap H)|<N+q$, where $N$ is the number of points of $C$ where $v_{n+1}=0$ and $q$ is the number of unbounded components of $C$.

Proof.

(of Theorem 8)  By induction on $k$. If $k=0$, there are no exponentials and the bound (12) reduces to the product of the degrees $d_{1}\cdots d_{n}$. This is the well known Bézout bound for a multivariate system of polynomial equations.

Now we will give the sketch of the proof and mention how our estimates change if the smooth assumptions in the induction step do not hold. In any case, the final inequalities needed to prove the bound (12) will hold.

For $k\geq 1$, to reduce the number of exponentials to $k-1$, we introduce a new variable $t$ such that the system with equations

[TABLE]

has the same as the original system when intersecting with the hyperplane $t=1$. We assume the functions (13) have [math] as a regular value so that the locus $\widehat{g_{1}}=\ldots=\widehat{g_{n}}=0$ is a smooth curve $C$ in $\mathbb{R}^{n+1}$ (this is a critical step that needs to be modified later), so we can apply the Khovanskii-Rolle Lemma. Indeed, by Cramer’s rule, the vector field

[TABLE]

is orthogonal to $\nabla\hat{g}_{i}$ for all $i$. So it is tangential to $C$ and non-vanishing because [math] is a regular value. Thus, the bound for $N$ is the number of solutions to the system in the $n+1$ variables $x_{1},x_{2},\ldots,x_{n},u$ with $u=ty_{k}(x)$ and $k-1$ exponentials

[TABLE]

Now $\frac{\partial y_{j}}{\partial x_{i}}(x)=y_{j}\cdot l_{ij}(x)$ where $l_{ij}$ is a linear function. Hence, $\frac{\partial\widehat{g_{i}}}{\partial x_{j}}=h_{ij}(x,y_{1}(x),\ldots,y_{k-1}(x),u)$ where $h_{ij}$ is a polynomial of degree at most $d_{i}+1$. Thus $v_{n+1}(x,u)$ is a polynomial of degree at most $(d_{1}+1)+\ldots+(d_{n}+1)=n+D$, where $D=d_{1}+\ldots+d_{n}$. By induction hypothesis,

[TABLE]

In order to bound $q$, the number of unbounded components of $C$, one observes that a hyperplane sufficiently far from the origin will meet $C$ in at least $q$ points (cf. [4, Lemma 12.6]). In other words, $q$ can be bounded by the number of solutions of a system

[TABLE]

for some $\lambda_{i},\mu\in\mathbb{R}$. Under non-degeneracy of the solutions, we get by induction hypothesis,

[TABLE]

So, in total,

[TABLE]

as we wanted. This is the end of the sketch.

If the smoothness assumptions for the system after introducing $t$ are not satisfied, the argument is modified via the Morse-Sard Theorem. The details of such modifications can be found along [11], although we find that for our theorem these are better summarized in [4, p. 293-295]. Essentially, one slightly perturbs the system from $\hat{g_{i}}=0$ to $\hat{g_{i}}=\epsilon_{i}$ ($\epsilon$ in a neighborhood of [math]) to guarantee obtaining a smooth curve $C$. The asserted bound (12) remains unchanged, and the number of non-degenerate solutions of the perturbed system cannot be less than the number for the original system.

Another change is that since the polynomial system might not define a proper map, one adds an extra variable $x_{0}$ with an extra equation

[TABLE]

with $R>0$ so that every preimage is now bounded. Morse-Sard now applies to conclude the set of regular values of $g=(g_{0},\ldots,g_{n}):\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n+1}$ is open and dense. The number of non-degenerate solutions of the new system has twice the number of non-degenerate solutions of the original system that lie in the open ball of radius $R$ centered at the origin (because from (19), a solution $(x_{1},x_{2},\ldots,x_{n})$ gives two possible values for $x_{0}$). In terms of bounding the corresponding $N^{\prime},q^{\prime}$, we now have

[TABLE]

(the extra 2 comes from the degree of (19), and we now have $n+2$ variables). For $q^{\prime}$, the bound becomes

[TABLE]

(the extra 1 from the hyperplane equation). Since the $R>0$ does not affect the bound computation, it can be taken large enough to include all the solutions to the original system. Thus $N^{\prime}+q^{\prime}$ is a bound for twice as many non-degenerate solutions of said original system. Finally, this way the induction step inequality can again be completed

[TABLE]

as needed. ∎

Now we can obtain Theorem 7 as a corollary of the above.

Proof.

(of Theorem 7) Let $\Phi(x)=\sum_{i=1}^{k}\alpha_{i}\Phi_{i}(x)$ be a Gaussian mixture and consider the system $g_{i}(x_{1},\ldots,x_{d})$ given by the partial derivatives $g_{i}=\frac{\partial\Phi}{\partial x_{i}}$. These can be interpreted as polynomials $F_{i}(x_{1},\ldots,x_{d},y_{1},\ldots,y_{k})$ by taking $y_{i}$ as the exponential kernel of $\Phi_{i}$. The system now has the form as in Theorem 8, and note that the degree of each $F_{i}$ is $2$. The number of non-degenerate critical points of $\Phi$ is thus the number of non-degenerate solutions to the system of $g_{i}$, and according to (12), it is bounded by

[TABLE]

∎

Finally, we show as promised that (11) is an upper bound on the number of modes of a Gaussian mixture, provided it is finite.

Theorem 10**.**

If a mixture of $k$ Gaussians in $\mathbb{R}^{d}$ has finitely many modes, then their number is bounded by

[TABLE]

Proof.

Let $\Phi(x)=\sum_{i=1}^{k}\alpha_{i}\Phi_{i}(x)$ be the pdf of the Gaussian mixture. If all of its modes are non-degenerate, then Theorem 8 applies and we’re done. The difficulty stems from considering possible degenerate modes. Note that by the finiteness hypothesis, all of them are isolated so we may fix disjoint neighborhoods $Q_{i}\subset\mathbb{R}^{d}$ over which each mode is the unique global maximum.

For any linear function $\ell(x_{1},\ldots,x_{d})=c\cdot x$, the function $\Phi+\ell$ has a gradient that differs by the constant vector $c$ from $\nabla\Phi$. In particular, the system given by the partial derivatives of $\Phi+\ell$ are still polynomials $F_{i}(x_{1},\ldots,x_{d},y_{1},\ldots,y_{k})$ of degree 2. By Theorem 7 we have the bound (11) on the non-degenerate modes of $\Phi+\ell$. Since $\Phi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is smooth, one of Morse’s Lemmas [12, Lemma A, p.11] states that for almost all $c\in\mathbb{R}^{d}$ (all except for a set of measure zero), $\Phi+\ell$ has only non-degenerate critical points.

Now, take $c\in\mathbb{R}^{d}$ in the complement of such measure zero set, with norm small enough so that $\Phi+\ell$ still has modes inside each of the neighborhoods $Q_{i}$. Then, $\Phi+\ell$ may have fewer critical points than $\Phi$ but we know $\Phi+\ell$ has at least as many modes as $\Phi$ does. Since the modes of $\Phi+\ell$ are bounded by (11), the result follows. ∎

Corollary 11**.**

If every mixture of $k$ Gaussians in $\mathbb{R}^{d}$ has finitely many modes, then

[TABLE]

6. Conclusion and Future Work

One of the motivations to study critical points of Gaussian mixtures comes from the mean shift algorithm. It converges if there are only finitely many critical points. This was the main goal sought in [17] for Gaussian kernels. As our bound (12) only bounds the number of non-degenerate critical points of Gaussian mixtures, a final answer is still open. However, since non-degenerate critical points are isolated, we see that the set of critical points of a Gaussian mixture is finite if and only if it consists of isolated points, since this set is closed and bounded (compare [8]).

Quantitatively, we do not expect our upper bound to be tight. Rather, proving the lower bound $\binom{d+k-1}{d}$ for all $d,k$ will be the main focus of a forthcoming paper, extending the technique used to prove Theorem 6.

We observe that the construction strategies used in our lower bound can be extended to elliptical distributions, not only Gaussians, and this could be pursued further. For example, an extension of Ray and Lindsay’s concept of ridgeline for mixtures of elliptical distributions is done in [1], including a study of modes for mixtures of two $t$-distributions.

Acknowledgements

The authors would like to thank Bernd Sturmfels for bringing the problem to their attention and for his supportive comments. Thanks to Peter Green for pointing to the discussion thread in the ANZstat mailing list. Carlos Améndola was supported by the Einstein Foundation Berlin. Alexander Engström would like to thank Günter Ziegler and Freie Universität Berlin for their hospitality during his sabbatical. The authors also express their gratitude to anonymous referees for their helpful suggestions.