Guaranteed Bounds on Information-Theoretic Measures of Univariate Mixtures Using Piecewise Log-Sum-Exp Inequalities

Abstract

Information-theoretic measures, such as the entropy, the cross-entropy and the Kullback–Leibler divergence between two mixture models, are core primitives in many signal processing tasks. Since the Kullback–Leibler divergence of mixtures provably does not admit a closed-form formula, it is in practice either estimated using costly Monte Carlo stochastic integration, approximated or bounded using various techniques. We present a fast and generic method that builds algorithmically closed-form lower and upper bounds on the entropy, the cross-entropy, the Kullback–Leibler and the α-divergences of mixtures. We illustrate the versatile method by reporting our experiments for approximating the Kullback–Leibler and the α-divergences between univariate exponential mixtures, Gaussian mixtures, Rayleigh mixtures and Gamma mixtures.

1. Introduction

Mixture models are commonly used in signal processing. A typical scenario is to use mixture models [1,2,3] to smoothly model histograms. For example, Gaussian Mixture Models (GMMs) can be used to convert grey-valued images into binary images by building a GMM fitting the image intensity histogram and then choosing the binarization threshold as the average of the Gaussian means [1]. Similarly, Rayleigh Mixture Models (RMMs) are often used in ultrasound imagery [2] to model histograms, and perform segmentation by classification. When using mixtures, a fundamental primitive is to define a proper statistical distance between them. The Kullback–Leibler (KL) divergence [4], also called relative entropy or information discrimination, is the most commonly used distance. Hence the main target of this paper is to faithfully measure the KL divergence. Let $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_i\left(x\right)$ and $m^{\prime }\left(x\right)={\sum }_{i=1}^{k^{\prime }}{w}_i^{\prime }{p}_i^{\prime }\left(x\right)$ be two finite statistical density mixtures of k and $k^{\prime }$ components, respectively. Notice that the Cumulative Density Function (CDF) of a mixture is like its density also a convex combinations of the component CDFs. However, beware that a mixture is not a sum of random variables (RVs). Indeed, sums of RVs have convolutional densities. In statistics, the mixture components $p_i\left(x\right)$ are often parametric: $p_i\left(x\right)=p(x;{\theta }_i)$, where ${\theta }_i$ is a vector of parameters. For example, a mixture of Gaussians (MoG also used as a shortcut instead of GMM) has each component distribution parameterized by its mean ${\mu }_i$ and its covariance matrix ${\mathsf{\Sigma }}_i$ (so that the parameter vector is ${\theta }_i=({\mu }_i,{\mathsf{\Sigma }}_i)$). Let $\mathcal{X}=\{x\in \mathbb{R}:p(x;\theta )>0\}$ be the support of the component distributions. Denote by $H_{\times }(m,m^{\prime })=-{\int }_{\mathcal{X}}m\left(x\right)log{m}^{\prime }\left(x\right)\mathrm{d}x$ the cross-entropy [4] between two continuous mixtures of densities m and $m^{\prime }$, and denote by $H\left(m\right)=H_{\times }(m,m)={\int }_{\mathcal{X}}m\left(x\right)log\frac{1}{m\left(x\right)}\mathrm{d}x=-{\int }_{\mathcal{X}}m\left(x\right)logm\left(x\right)\mathrm{d}x$ the Shannon entropy [4]. Then the Kullback–Leibler divergence between m and $m^{\prime }$ is given by:

$$\mathrm{KL}(m:m^{\prime })=H_{\times }(m,m^{\prime })-H\left(m\right)={\int }_{\mathcal{X}}m\left(x\right)log\frac{m\left(x\right)}{m^{\prime }\left(x\right)}\mathrm{d}x\ge 0.$$

(1)

The notation “:” is used instead of the usual comma “,” notation to emphasize that the distance is not a metric distance since neither is it symmetric ($\mathrm{KL}(m:m^{\prime })\ne \mathrm{KL}(m^{\prime }:m)$), nor does it satisfy the triangular inequality [4] of metric distances ($\mathrm{KL}(m:m^{\prime })+\mathrm{KL}(m^{\prime }:m^{″})\ngeqq \mathrm{KL}(m:m^{\prime \prime })$). When the natural base of the logarithm is chosen, we get a differential entropy measure expressed in nat units. Alternatively, we can also use the base-2 logarithm (${log}_{2}x=\frac{logx}{log2}$) and get the entropy expressed in bit units. Although the KL divergence is available in closed-form for many distributions (in particular as equivalent Bregman divergences for exponential families [5], see Appendix C), it was proven that the Kullback–Leibler divergence between two (univariate) GMMs is not analytic [6] (see also the particular case of a GMM of two components with the same variance that was analyzed in [7]). See Appendix A for an analysis. Note that the differential entropy may be negative. For example, the differential entropy of a univariate Gaussian distribution is $log\left(\sigma \sqrt{2\pi e}\right)$, and is therefore negative when the standard variance $\sigma <\frac{1}{\sqrt{2\pi e}}\approx 0.242$. We consider continuous distributions with entropies well-defined (entropy may be undefined for singular distributions like Cantor’s distribution [8]).

1.1. Prior Work

Many approximation techniques have been designed to beat the computationally intensive Monte Carlo (MC) stochastic estimation: ${\widehat{\mathrm{KL}}}_s(m:m^{\prime })=\frac{1}{s}{\sum }_{i=1}^{s}log\frac{m\left(x_i\right)}{m^{\prime }\left(x_i\right)}$ with $x_1,\dots ,x_s\sim m\left(x\right)$ (s independently and identically distributed (i.i.d.) samples $x_1,\dots ,x_s$). The MC estimator is asymptotically consistent, ${lim}_{s\to \infty }{\widehat{\mathrm{KL}}}_s(m:m^{\prime })=\mathrm{KL}(m:m^{\prime })$, so that the “true value” of the KL of mixtures is estimated in practice by taking a very large sample (say, $s={10}^9$). However, we point out that the MC estimator gives as output a stochastic approximation, and therefore does not guarantee deterministic bounds (confidence intervals may be used). Deterministic lower and upper bounds of the integral can be obtained by various numerical integration techniques using quadrature rules. We refer to [9,10,11,12] for the current state-of-the-art approximation techniques and bounds on the KL of GMMs. The latest work for computing the entropy of GMMs is [13]. It considers arbitrary finely tuned bounds of the entropy of isotropic Gaussian mixtures (a case encountered when dealing with KDEs, kernel density estimators). However, there is a catch in the technique of [13]: It relies on solving the unique roots of some log-sum-exp equations (See Theorem 1 of [13], p. 3342) that do not admit a closed-form solution. Thus it is a hybrid method that contrasts with our combinatorial approach. Bounds of the KL divergence between mixture models can be generalized to bounds of the likelihood function of mixture models [14], because log-likelihood is just the KL between the empirical distribution and the mixture model up to a constant shift.

In information geometry [15], a mixture family of linearly independent probability distributions $p_1\left(x\right),...,p_k\left(x\right)$ is defined by the convex combination of those non-parametric component distributions: $m(x;\eta )={\sum }_{i=1}^k{\eta }_i{p}_i\left(x\right)$ with ${\eta }_i>0$ and ${\sum }_{i=1}^k{\eta }_i=1$. A mixture family induces a dually flat space where the Kullback–Leibler divergence is equivalent to a Bregman divergence [5,15] defined on the η-parameters. However, in that case, the Bregman convex generator $F\left(\eta \right)=\int m(x;\eta )logm(x;\eta )\mathrm{d}x$ (the Shannon information) is not available in closed-form. Except for the family of multinomial distributions that is both a mixture family (with closed-form $\mathrm{KL}(m:m^{\prime })={\sum }_{i=1}^k{m}_{i}log\frac{m_i}{m_i^{\prime }}$, the discrete KL [4]) and an exponential family [15].

1.2. Contributions

In this work, we present a simple and efficient method that builds algorithmically a closed-form formula that guarantees both deterministic lower and upper bounds on the KL divergence within an additive factor of $logk+log{k}^{\prime }$. We then further refine our technique to get improved adaptive bounds. For univariate GMMs, we get the non-adaptive bounds in $O(klogk+k^{\prime }log{k}^{\prime })$ time, and the adaptive bounds in $O(k^2+k^{\prime 2})$ time. To illustrate our generic technique, we demonstrate it based on Exponential Mixture Models (EMMs), Gamma mixtures, RMMs and GMMs. We extend our preliminary results on KL divergence [16] to other information theoretical measures such as the differential entropy and α-divergences.

1.3. Paper Outline

The paper is organized as follows. Section 2 describes the algorithmic construction of the formula using piecewise log-sum-exp inequalities for the cross-entropy and the Kullback–Leibler divergence. Section 3 instantiates this algorithmic principle to the entropy and discusses related works. Section 4 extends the proposed bounds to the family of alpha divergences. Section 5 discusses an extension of the lower bound to f-divergences. Section 6 reports our experimental results on several mixture families. Finally, Section 7 concludes this work by discussing extensions to other statistical distances. Appendix A proves that the Kullback–Leibler divergence of mixture models is not analytic [6]. Appendix B reports the closed-form formula for the KL divergence between scaled and truncated distributions of the same exponential family [17] (that include Rayleigh, Gaussian and Gamma distributions among others). Appendix C shows that the KL divergence between two mixtures can be approximated by a Bregman divergence.

2. A Generic Combinatorial Bounding Algorithm Based on Density Envelopes

Let us bound the cross-entropy $H_{\times }(m:m^{\prime })$ by deterministic lower and upper bounds, $L_{\times }(m:m^{\prime })\le H_{\times }(m:m^{\prime })\le U_{\times }(m:m^{\prime })$, so that the bounds on the Kullback–Leibler divergence $\mathrm{KL}(m:m^{\prime })=H_{\times }(m:m^{\prime })-H_{\times }(m:m)$ follows as:

$$L_{\times }(m:m^{\prime })-U_{\times }(m:m)\le \mathrm{KL}(m:m^{\prime })\le U_{\times }(m:m^{\prime })-L_{\times }(m:m).$$

(2)

Since the cross-entropy of two mixtures ${\sum }_{i=1}^k{w}_i{p}_i\left(x\right)$ and ${\sum }_{j=1}^{k^{\prime }}{w}_j^{\prime }{p}_j^{\prime }\left(x\right)$:

$$H_{\times }(m:m^{\prime })=-{\int }_{\mathcal{X}}\left(\sum _{i=1}^k{w}_i{p}_i\left(x\right)\right)log\left(\sum _{j=1}^{k^{\prime }}{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\right)\mathrm{d}x$$

(3)

has a log-sum term of positive arguments, we shall use bounds on the log-sum-exp (lse) function [18,19]:

$$\mathrm{lse}\left({\left\{x_i\right\}}_{i=1}^l\right)=log\left(\sum _{i=1}^l{e}^{x_i}\right).$$

We have the following basic inequalities:

$$max{\left\{x_i\right\}}_{i=1}^l<\mathrm{lse}\left({\left\{x_i\right\}}_{i=1}^l\right)\le logl+max{\left\{x_i\right\}}_{i=1}^l.$$

(4)

The left-hand-side (LHS) strict inequality holds because ${\sum }_{i=1}^l{e}^{x_i}>max{\left\{e^{x_i}\right\}}_{i=1}^l=exp\left(max{\left\{x_i\right\}}_{i=1}^l\right)$ since $e^x>0,\forall x\in \mathbb{R}$. The right-hand-side (RHS) inequality follows from the fact that ${\sum }_{i=1}^l{e}^{x_i}\le lmax{\left\{e^{x_i}\right\}}_{i=1}^l=lexp(max{\left\{x_i\right\}}_{i=1}^l)$, and equality holds if and only if $x_1=\cdots =x_l$. The lse function is convex but not strictly convex, see exercise 7.9 [20]. It is known [21] that the conjugate of the lse function is the negative entropy restricted to the probability simplex. The lse function enjoys the following translation identity property: $\mathrm{lse}\left({\left\{x_i\right\}}_{i=1}^l\right)=c+\mathrm{lse}\left({\{x_i-c\}}_{i=1}^l\right)$, $\forall c\in \mathbb{R}$. Similarly, we can also lower bound the lse function by $logl+min{\left\{x_i\right\}}_{i=1}^l$. We write equivalently that for l positive numbers $x_1,\dots ,x_l$,

$$max\left\{logmax{\left\{x_i\right\}}_{i=1}^l,logl+logmin{\left\{x_i\right\}}_{i=1}^l\right\}\le log\sum _{i=1}^l{x}_i\le logl+logmax{\left\{x_i\right\}}_{i=1}^l.$$

(5)

In practice, we seek matching lower and upper bounds that minimize the bound gap. The gap of that ham-sandwich inequality in Equation (5) is $min\{log\frac{{max}_i{x}_i}{{min}_i{x}_i},logl\}$, which is upper bounded by $logl$.

A mixture model ${\sum }_{j=1}^{k^{\prime }}{w}_j^{\prime }{p}_j^{\prime }\left(x\right)$ must satisfy

$$\begin{array}{c}max\left\{max{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }},log{k}^{\prime }+min{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}\right\}\hfill \\ \le log\left(\sum _{j=1}^{k^{\prime }}{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\right)\le log{k}^{\prime }+max{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}\hfill \end{array}$$

(6)

point-wisely for any $x\in \mathcal{X}$. Therefore we shall bound the integral term ${\int }_{\mathcal{X}}m\left(x\right)log\left({\sum }_{j=1}^{k^{\prime }}{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\right)\mathrm{d}x$ in Equation (3) using piecewise lse inequalities where the min and max are kept unchanged. We get

$$\begin{array}{cc}\hfill L_{\times }(m:m^{\prime })& =-{\int }_{\mathcal{X}}m\left(x\right)max{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}\mathrm{d}x-log{k}^{\prime },\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill U_{\times }(m:m^{\prime })& =-{\int }_{\mathcal{X}}m\left(x\right)max\left\{min{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}+log{k}^{\prime },max{\{log{w}_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}\right\}\mathrm{d}x.\hfill \end{array}$$

(8)

In order to calculate $L_{\times }(m:m^{\prime })$ and $U_{\times }(m:m^{\prime })$ efficiently using closed-form formula, let us compute the upper and lower envelopes of the $k^{\prime }$ real-valued functions ${\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$ defined on the support $\mathcal{X}$, that is, ${\mathcal{E}}_{\mathrm{U}}\left(x\right)=max{\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$ and ${\mathcal{E}}_{\mathrm{L}}\left(x\right)=min{\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$. These envelopes can be computed exactly using techniques of computational geometry [22,23] provided that we can calculate the roots of the equation $w_r^{\prime }{p}_r^{\prime }\left(x\right)=w_s^{\prime }{p}_s^{\prime }\left(x\right)$, where $w_r^{\prime }{p}_r^{\prime }\left(x\right)$ and $w_s^{\prime }{p}_s^{\prime }\left(x\right)$ are a pair of weighted components. (Although this amounts to solve quadratic equations for Gaussian or Rayleigh distributions, the roots may not always be available in closed form, e.g. in the case of Weibull distributions.)

Let the envelopes be combinatorially described by ℓ elementary interval pieces in the form $I_r=(a_r,a_{r+1})$ partitioning the support $\mathcal{X}={\uplus }_{r=1}^{\ell }{I}_r$ (with $a_1=min\mathcal{X}$ and $a_{\ell +1}=max\mathcal{X}$). Observe that on each interval $I_r$, the maximum of the functions ${\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$ is given by $w_{\delta \left(r\right)}^{\prime }{p}_{\delta \left(r\right)}^{\prime }\left(x\right)$, where $\delta \left(r\right)$ indicates the weighted component dominating all the others, i.e., the arg max of ${\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$ for any $x\in I_r$, and the minimum of ${\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}_{j=1}^{k^{\prime }}$ is given by $w_{ϵ\left(r\right)}^{\prime }{p}_{ϵ\left(r\right)}^{\prime }\left(x\right)$.

To fix ideas, when mixture components are univariate Gaussians, the upper envelope ${\mathcal{E}}_{\mathrm{U}}\left(x\right)$ amounts to find equivalently the lower envelope of $k^{\prime }$ parabolas (see Figure 1) which has linear complexity, and can be computed in $O(k^{\prime }log{k}^{\prime })$-time [24], or in output-sensitive time $O(k^{\prime }log\ell )$ [25], where ℓ denotes the number of parabola segments in the envelope. When the Gaussian mixture components have all the same weight and variance (e.g., kernel density estimators), the upper envelope amounts to find a lower envelope of cones: ${min}_j|x-{\mu }_j^{\prime }|$ (a Voronoi diagram in arbitrary dimension).

To proceed once the envelopes have been built, we need to calculate two types of definite integrals on those elementary intervals: (i) the probability mass in an interval ${\int }_a^{b}p\left(x\right)\mathrm{d}x=\mathrm{\Phi }\left(b\right)-\mathrm{\Phi }\left(a\right)$ where Φ denotes the Cumulative Distribution Function (CDF); and (ii) the partial cross-entropy $-{\int }_a^{b}p\left(x\right)log{p}^{\prime }\left(x\right)\mathrm{d}x$ [26]. Thus let us define these two quantities:

$$\begin{array}{ccc}\hfill C_{i,j}(a,b)& =& -{\int }_a^b{w}_i{p}_i\left(x\right)log(w_j^{\prime }{p}_j^{\prime }\left(x\right))\mathrm{d}x,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill M_i(a,b)& =& -{\int }_a^b{w}_i{p}_i\left(x\right)\mathrm{d}x.\hfill \end{array}$$

(10)

By Equations (7) and (8), we get the bounds of $H_{\times }(m:m^{\prime })$ as

$$\begin{array}{cc}\hfill L_{\times }(m:m^{\prime })& =\sum _{r=1}^{\ell }\sum _{s=1}^k{C}_{s,\delta \left(r\right)}(a_r,a_{r+1})-log{k}^{\prime },\hfill \\ \hfill U_{\times }(m:m^{\prime })& =\sum _{r=1}^{\ell }\sum _{s=1}^{k}min\left\{C_{s,\delta \left(r\right)}(a_r,a_{r+1}),C_{s,ϵ\left(r\right)}(a_r,a_{r+1})-M_s(a_r,a_{r+1})log{k}^{\prime }\right\}.\hfill \end{array}$$

(11)

The size of the lower/upper bound formula depends on the envelope complexity ℓ, the number k of mixture components, and the closed-form expressions of the integral terms $C_{i,j}(a,b)$ and $M_i(a,b)$. In general, when a pair of weighted component densities intersect in at most p points, the envelope complexity is related to the Davenport–Schinzel sequences [27]. It is quasi-linear for bounded $p=O\left(1\right)$, see [27].

Note that in symbolic computing, the Risch semi-algorithm [28] solves the problem of computing indefinite integration in terms of elementary functions provided that there exists an oracle (hence the term “semi-algorithm”) for checking whether an expression is equivalent to zero or not (however it is unknown whether there exists an algorithm implementing the oracle or not).

We presented the technique by bounding the cross-entropy (and entropy) to deliver lower/upper bounds on the KL divergence. When only the KL divergence needs to be bounded, we rather consider the ratio term $\frac{m\left(x\right)}{m^{\prime }\left(x\right)}$. This requires to partition the support $\mathcal{X}$ into elementary intervals by overlaying the critical points of both the lower and upper envelopes of $m\left(x\right)$ and $m^{\prime }\left(x\right)$, which can be done in linear time. In a given elementary interval, since $max\{k{min}_i\left\{w_i{p}_i\left(x\right)\right\},{max}_i\left\{w_i{p}_i\left(x\right)\right\}\}\le m\left(x\right)\le k{max}_i\left\{w_i{p}_i\left(x\right)\right\}$, we then consider the inequalities:

$$\begin{array}{c}\hfill \frac{max\{k{min}_i\left\{w_i{p}_i\left(x\right)\right\},{max}_i\left\{w_i{p}_i\left(x\right)\right\}\}}{k^{\prime }{max}_j\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}}\le \frac{m\left(x\right)}{m^{\prime }\left(x\right)}\le \frac{k{max}_i\left\{w_i{p}_i\left(x\right)\right\}}{max\{k^{\prime }{min}_j\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\},{max}_j\{w_j^{\prime }{p}_j^{\prime }\left(x\right)\}\}}.\end{array}$$

(12)

We now need to compute definite integrals of the form ${\int }_a^b{w}_{1}p(x;{\theta }_1)log\frac{w_{2}p(x;{\theta }_2)}{w_{3}p(x;{\theta }_3)}\mathrm{d}x$ (see Appendix B for explicit formulas when considering scaled and truncated exponential families [17]). (Thus for exponential families, the ratio of densities removes the auxiliary carrier measure term.)

We call these bounds CELB and CEUB for Combinatorial Envelope Lower and Upper Bounds, respectively.

2.1. Tighter Adaptive Bounds

We shall now consider shape-dependent bounds improving over the additive $logk+log{k}^{\prime }$ non-adaptive bounds. This is made possible by a decomposition of the lse function explained as follows. Let $t_i(x_1,\dots ,x_k)=log\left({\sum }_{j=1}^k{e}^{x_j-x_i}\right)$. By translation identity of the lse function,

$$\mathrm{lse}(x_1,\dots ,x_k)=x_i+t_i(x_1,\dots ,x_k)$$

(13)

for all $i\in \left[k\right]$. Since $e^{x_j-x_i}=1$ if $j=i$, and $e^{x_j-x_i}>0$, we have necessarily $t_i(x_1,\dots ,x_k)>0$ for any $i\in \left[k\right]$. Since Equation (13) is an identity for all $i\in \left[k\right]$, we minimize the residual $t_i(x_1,\dots ,x_k)$ by maximizing $x_i$. Denoting by $x_{\left(1\right)},\dots ,x_{\left(k\right)}$ the sequence of numbers sorted in non-decreasing order, the decomposition

$$\mathrm{lse}(x_1,\dots ,x_k)=x_{\left(k\right)}+t_{\left(k\right)}(x_1,\dots ,x_k)$$

(14)

yields the smallest residual. Since $x_{\left(j\right)}-x_{\left(k\right)}\le 0$ for all $j\in \left[k\right]$, we have

$$t_{\left(k\right)}\left(x_1,\dots ,x_k\right)=log\left(1+\sum _{j=1}^{k-1}{e}^{x_{\left(j\right)}-x_{\left(k\right)}}\right)\le logk.$$

This shows the bounds introduced earlier can indeed be improved by a more accurate computation of the residual term $t_{\left(k\right)}\left(x_1,\dots ,x_k\right)$.

When considering 1D GMMs, let us now bound $t_{\left(k\right)}(x_1,\dots ,x_k)$ in a combinatorial range $I_r=(a_r,a_{r+1})$. Let $\delta =\delta \left(r\right)$ denote the index of the dominating weighted component in this range. Then,

$$\forall x\in I_r,\forall i,exp\left(-log{\sigma }_i-\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}+log{w}_i\right)\le exp\left(-log{\sigma }_{\delta }-\frac{{(x-{\mu }_{\delta })}^2}{2{\sigma }_{\delta }^2}+log{w}_{\delta }\right).$$

Thus we have:

$$\begin{array}{c}\hfill logm\left(x\right)=log\frac{w_{\delta }}{{\sigma }_{\delta }\sqrt{2\pi }}-\frac{{(x-{\mu }_{\delta })}^2}{2{\sigma }_{\delta }^2}+log\left(1+\sum _{i\ne \delta }exp\left(-\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}+log\frac{w_i}{{\sigma }_i}+\frac{{(x-{\mu }_{\delta })}^2}{2{\sigma }_{\delta }^2}-log\frac{w_{\delta }}{{\sigma }_{\delta }}\right)\right).\end{array}$$

Now consider the ratio term:

$${\rho }_{i,\delta }\left(x\right)=exp\left(-\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}+log\frac{w_i{\sigma }_{\delta }}{w_{\delta }{\sigma }_i}+\frac{{(x-{\mu }_{\delta })}^2}{2{\sigma }_{\delta }^2}\right).$$

It is maximized in $I_r=(a_r,a_{r+1})$ by maximizing equivalently the following quadratic equation:

$$l_{i,\delta }\left(x\right)=-\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}+log\frac{w_i{\sigma }_{\delta }}{w_{\delta }{\sigma }_i}+\frac{{(x-{\mu }_{\delta })}^2}{2{\sigma }_{\delta }^2}.$$

Setting the derivative to zero ($l_{i,\delta }^{\prime }\left(x\right)=0$), we get the root (when ${\sigma }_i\ne {\sigma }_{\delta }$)

$$x_{i,\delta }=\left(\frac{{\mu }_{\delta }}{{\sigma }_{\delta }^2}-\frac{{\mu }_i}{{\sigma }_i^2}\right)/\left(\frac{1}{{\sigma }_{\delta }^2}-\frac{1}{{\sigma }_i^2}\right).$$

If $x_{i,\delta }\in I_r$, the ratio ${\rho }_{i,\delta }\left(x\right)$ can be bounded in the slab $I_r$ by considering the extreme values of the three element set $\{{\rho }_{i,\delta }\left(a_r\right),{\rho }_{i,\delta }\left(x_{i,\delta }\right),{\rho }_{i,\delta }\left(a_{r+1}\right)\}$. Otherwise ${\rho }_{i,\delta }\left(x\right)$ is monotonic in $I_r$, its bounds in $I_r$ are given by $\{{\rho }_{i,\delta }\left(a_r\right),{\rho }_{i,\delta }\left(a_{r+1}\right)\}$. In any case, let ${\rho }_{i,\delta }^{min}\left(r\right)$ and ${\rho }_{i,\delta }^{max}\left(r\right)$ represent the resulting lower and upper bounds of ${\rho }_{i,\delta }\left(x\right)$ in $I_r$. Then $t_{\delta }$ is bounded in the range $I_r$ by:

$$0<log\left(1+\sum _{i\ne \delta }{\rho }_{i,\delta }^{min}\left(r\right)\right)\le t_{\delta }\le log\left(1+\sum _{i\ne \delta }{\rho }_{i,\delta }^{max}\left(r\right)\right)\le logk.$$

In practice, we always get better bounds using the shape-dependent technique at the expense of computing overall $O\left(k^2\right)$ intersection points of the pairwise densities. We call those bounds CEALB and CEAUB for Combinatorial Envelope Adaptive Lower Bound and Combinatorial Envelope Adaptive Upper Bound.

Let us illustrate one scenario where this adaptive technique yields very good approximations. Consider a GMM with all variance ${\sigma }^2$ tending to zero (a mixture of k Diracs). Then in a combinatorial slab $I_r$, we have ${\rho }_{i,\delta }^{max}\left(r\right)\to 0$ for all $i\ne \delta$, and therefore we get tight bounds.

As a related technique, we could also upper bound ${\int }_{a_r}^{a_{r+1}}logm\left(x\right)\mathrm{d}x$ by $(a_{r+1}-a_r)logm(a_r,a_{r+1})$ where $m(x,x^{\prime })$ denotes the maximal value of the mixture density in the range $(x,x^{\prime })$. This maximal value is either found at the slab extremities, or is a mode of the GMM. It then requires to find the modes of a GMM [29,30], for which no analytical solution is known in general.

2.2. Another Derivation Using the Arithmetic-Geometric Mean Inequality

Let us start by considering the inequality of arithmetic and geometric weighted means (AGI, Arithmetic-Geometric Inequality) applied to the mixture component distributions:

$$m\left(x\right)=\sum _{i=1}^k{w}_{i}p(x;{\theta }_i)\ge \prod _{i=1}^{k}p{(x;{\theta }_i)}^{w_i}$$

with equality holds iff. ${\theta }_1=\dots ={\theta }_k$.

To get a tractable formula with a positive remainder of the log-sum term $logm\left(x\right)$, we need to have the log argument greater or equal to 1, and thus we shall write the positive remainder:

$$R\left(x\right)=log\left(\frac{m\left(x\right)}{{\prod }_{i=1}^{k}p{(x;{\theta }_i)}^{w_i}}\right)\ge 0.$$

Therefore, we can decompose the log-sum into a tractable part and a remainder as:

$$logm\left(x\right)=\sum _{i=1}^k{w}_{i}logp(x;{\theta }_i)+log\left(\frac{m\left(x\right)}{{\prod }_{i=1}^{k}p{(x;{\theta }_i)}^{w_i}}\right).$$

(15)

For exponential families, the first term can be integrated accurately. For the second term, we notice that ${\prod }_{i=1}^{k}p{(x;{\theta }_i)}^{w_i}$ is a distribution in the same exponential family. We denote $p(x;{\theta }_0)={\prod }_{i=1}^{k}p{(x;{\theta }_i)}^{w_i}$. Then

$$R\left(x\right)=log\left(\sum _{i=1}^k{w}_i\frac{p(x;{\theta }_i)}{p(x;{\theta }_0)}\right)$$

As the ratio $p(x;{\theta }_i)/p(x;{\theta }_0)$ can be bounded above and below using techniques in Section 2.1, $R\left(x\right)$ can be correspondingly bounded. Notice the similarity between Equations (14) and (15). The key difference with the adaptive bounds is that, here we choose $p(x;{\theta }_0)$ instead of the dominating component in $m\left(x\right)$ as the “reference distribution” in the decomposition. This subtle difference is not presented in detail in our experimental studies but discussed here for completeness. Essentially, the gap of the bounds is up to the difference between the geometric average and the arithmetic average. In the extreme case that all mixture components are identical, this gap will reach zero. Therefore we expect good quality bounds with a small gap when the mixture components are similar as measured by KL divergence.

2.3. Case Studies

In the following, we instantiate the proposed method for several prominent cases on the mixture of exponential family distributions.

2.3.1. The Case of Exponential Mixture Models

An exponential distribution has density $p(x;\lambda )=\lambda exp(-\lambda x)$ defined on $\mathcal{X}=[0,\infty )$ for $\lambda >0$. Its CDF is $\mathrm{\Phi }(x;\lambda )=1-exp(-\lambda x)$. Any two components $w_{1}p(x;{\lambda }_1)$ and $w_{2}p(x;{\lambda }_2)$ (with ${\lambda }_1\ne {\lambda }_2$) have a unique intersection point

$$x^{\star }=\frac{log\left(w_1{\lambda }_1\right)-log\left(w_2{\lambda }_2\right)}{{\lambda }_1-{\lambda }_2}$$

(16)

if $x^{\star }\ge 0$; otherwise they do not intersect. The basic formulas to evaluate the bounds are

$$\begin{array}{cc}\hfill C_{i,j}(a,b)=& log\left({\lambda }_j^{\prime }{w}_j^{\prime }\right)M_i(a,b)+w_i{\lambda }_j^{\prime }\left[\left(a+\frac{1}{{\lambda }_i}\right)e^{-{\lambda }_{i}a}-\left(b+\frac{1}{{\lambda }_i}\right)e^{-{\lambda }_{i}b}\right],\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill M_i(a,b)=& -w_i\left(e^{-{\lambda }_{i}a}-e^{-{\lambda }_{i}b}\right).\hfill \end{array}$$

(18)

2.3.2. The Case of Rayleigh Mixture Models

A Rayleigh distribution has density $p(x;\sigma )=\frac{x}{{\sigma }^2}exp\left(-\frac{x^2}{2{\sigma }^2}\right)$, defined on $\mathcal{X}=[0,\infty )$ for $\sigma >0$. Its CDF is $\mathrm{\Phi }(x;\sigma )=1-exp\left(-\frac{x^2}{2{\sigma }^2}\right)$. Any two components $w_{1}p(x;{\sigma }_1)$ and $w_{2}p(x;{\sigma }_2)$ (with ${\sigma }_1\ne {\sigma }_2$) must intersect at $x_0=0$ and can have at most one other intersection point

$$x^{\star }=\sqrt{log\frac{w_1{\sigma }_2^2}{w_2{\sigma }_1^2}/\left(\frac{1}{2{\sigma }_1^2}-\frac{1}{2{\sigma }_2^2}\right)}$$

(19)

if the square root is well defined and $x^{\star }>0$. We have

$$\begin{array}{cc}\hfill C_{i,j}(a,b)=& log\frac{w_j^{\prime }}{{({\sigma }_j^{\prime })}^2}{M}_i(a,b)+\frac{w_i}{2{({\sigma }_j^{\prime })}^2}\left[(a^2+2{\sigma }_i^2)e^{-\frac{a^2}{2{\sigma }_i^2}}-(b^2+2{\sigma }_i^2)e^{-\frac{b^2}{2{\sigma }_i^2}}\right]\hfill \\ & -w_i{\int }_a^b\frac{x}{{\sigma }_i^2}exp\left(-\frac{x^2}{2{\sigma }_i^2}\right)logx\mathrm{d}x,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill M_i(a,b)=& -w_i\left(e^{-\frac{a^2}{2{\sigma }_i^2}}-e^{-\frac{b^2}{2{\sigma }_i^2}}\right).\hfill \end{array}$$

(21)

The last term in Equation (20) does not have a simple closed form (it requires the exponential integral, Ei). One need a numerical integrator to compute it.

2.3.3. The Case of Gaussian Mixture Models

The Gaussian density $p(x;\mu ,\sigma )=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-{(x-\mu )}^2/\left(2{\sigma }^2\right)}$ has support $\mathcal{X}=\mathbb{R}$ and parameters $\mu \in \mathbb{R}$ and $\sigma >0$. Its CDF is $\mathrm{\Phi }(x;\mu ,\sigma )=\frac{1}{2}\left[1+\mathrm{erf}(\frac{x-\mu }{\sqrt{2}\sigma })\right]$, where $\mathrm{erf}$ is the Gauss error function. The intersection point $x^{\star }$ of two components $w_{1}p(x;{\mu }_1,{\sigma }_1)$ and $w_{2}p(x;{\mu }_2,{\sigma }_2)$ can be obtained by solving the quadratic equation $log\left(w_{1}p(x;{\mu }_1,{\sigma }_1)\right)=log\left(w_{2}p(x;{\mu }_2,{\sigma }_2)\right)$, which gives at most two solutions. As shown in Figure 1, the upper envelope of Gaussian densities corresponds to the lower envelope of parabolas. We have

$$\begin{array}{cc}\hfill C_{i,j}(a,b)=& M_i(a,b)\left(log{w}_j^{\prime }-log{\sigma }_j^{\prime }-\frac{1}{2}log\left(2\pi \right)-\frac{1}{2{({\sigma }_j^{\prime })}^2}\left({({\mu }_j^{\prime }-{\mu }_i)}^2+{\sigma }_i^2\right)\right)\hfill \\ & +\frac{w_i{\sigma }_i}{2\sqrt{2\pi }{({\sigma }_j^{\prime })}^2}\left[(a+{\mu }_i-2{\mu }_j^{\prime })e^{-\frac{{(a-{\mu }_i)}^2}{2{\sigma }_i^2}}-(b+{\mu }_i-2{\mu }_j^{\prime })e^{-\frac{{(b-{\mu }_i)}^2}{2{\sigma }_i^2}}\right],\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill M_i(a,b)=& -\frac{w_i}{2}\left(\mathrm{erf}\left(\frac{b-{\mu }_i}{\sqrt{2}{\sigma }_i}\right)-\mathrm{erf}\left(\frac{a-{\mu }_i}{\sqrt{2}{\sigma }_i}\right)\right).\hfill \end{array}$$

(23)

2.3.4. The Case of Gamma Distributions

For simplicity, we only consider gamma distributions with the shape parameter $k>0$ fixed and the scale $\lambda >0$ varying. The density is defined on $(0,\infty )$ as $p(x;k,\lambda )=\frac{x^{k-1}{e}^{-\frac{x}{\lambda }}}{{\lambda }^k\mathrm{\Gamma }\left(k\right)}$, where $\mathrm{\Gamma }(·)$ is the gamma function. Its CDF is $\mathrm{\Phi }(x;k,\lambda )=\gamma (k,x/\lambda )/\mathrm{\Gamma }\left(k\right)$, where $\gamma (·,·)$ is the lower incomplete gamma function. Two weighted gamma densities $w_{1}p(x;k,{\lambda }_1)$ and $w_{2}p(x;k,{\lambda }_2)$ (with ${\lambda }_1\ne {\lambda }_2$) intersect at a unique point

$$x^{\star }=\left(log\frac{w_1}{{\lambda }_1^k}-log\frac{w_2}{{\lambda }_2^k}\right)/\left(\frac{1}{{\lambda }_1}-\frac{1}{{\lambda }_2}\right)$$

(24)

if $x^{\star }>0$; otherwise they do not intersect. From straightforward derivations,

$$\begin{array}{cc}\hfill C_{i,j}(a,b)& =log\frac{w_j^{\prime }}{{({\lambda }_j^{\prime })}^k\mathrm{\Gamma }\left(k\right)}{M}_i(a,b)+w_i{\int }_a^b\frac{x^{k-1}{e}^{-\frac{x}{{\lambda }_i}}}{{\lambda }_i^k\mathrm{\Gamma }\left(k\right)}\left(\frac{x}{{\lambda }_j^{\prime }}-(k-1)logx\right)dx,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill M_i(a,b)& =-\frac{w_i}{\mathrm{\Gamma }\left(k\right)}\left(\gamma \left(k,\frac{b}{{\lambda }_i}\right)-\gamma \left(k,\frac{a}{{\lambda }_i}\right)\right).\hfill \end{array}$$

(26)

Similar to the case of Rayleigh mixtures, the last term in Equation (25) relies on numerical integration.

3. Upper-Bounding the Differential Entropy of a Mixture

First, consider a finite parametric mixture model $m\left(x\right)={\sum }_{i=1}^k{w}_{i}p(x;{\theta }_i)$. Using the chain rule of the entropy, we end up with the well-known lemma:

Lemma 1.

The entropy of a d-variate mixture is upper bounded by the sum of the entropy of its marginal mixtures: $H\left(m\right)\le {\sum }_{i=1}^{d}H\left(m_i\right)$, where $m_i$ is the 1D marginal mixture with respect to variable $x_i$.

Since the 1D marginals of a multivariate GMM are univariate GMMs, we thus get a loose upper bound. A generic sample-based probabilistic bound is reported for the entropies of distributions with given support [31]: The method builds probabilistic upper and lower piecewisely linear CDFs based on an i.i.d. finite sample set of size n and a given deviation probability threshold. It then builds algorithmically between those two bounds the maximum entropy distribution [31] with a so-called string-tightening algorithm.

Instead, we proceed as follows: Consider finite mixtures of component distributions defined on the full support ${\mathbb{R}}^d$ that have finite component means and variances (like exponential families). Then we shall use the fact that the maximum entropy distribution with prescribed mean and variance is a Gaussian distribution, and conclude the upper bound by plugging the mixture mean and variance in the differential entropy formula of the Gaussian distribution. In general, the maximum entropy with moment constraints yields as a solution an exponential family.

Without loss of generality, consider GMMs in the form $m\left(x\right)={\sum }_{i=1}^k{w}_{i}p(x;{\mu }_i,{\mathsf{\Sigma }}_i)$ (${\mathsf{\Sigma }}_i={\sigma }_i^2$ for univariate Gaussians). The mean $\overline{\mu }$ of the mixture is $\overline{\mu }={\sum }_{i=1}^k{w}_i{\mu }_i$ and the variance is ${\overline{\sigma }}^2=E\left[m^2\right]-E{\left[m\right]}^2$. Since $E\left[m^2\right]={\sum }_{i=1}^k{w}_i\int x^{2}p(x;{\mu }_i,{\mathsf{\Sigma }}_i)\mathrm{d}x={\sum }_{i=1}^k{w}_i\left({\mu }_i^2+{\sigma }_i^2\right)$, we deduce that

$${\overline{\sigma }}^2=\sum _{i=1}^k{w}_i({\mu }_i^2+{\sigma }_i^2)-{\left(\sum _{i=1}^k{w}_i{\mu }_i\right)}^2=\sum _{i=1}^k{w}_i\left[{({\mu }_i-\overline{\mu })}^2+{\sigma }_i^2\right].$$

The entropy of a random variable with a prescribed variance ${\overline{\sigma }}^2$ is maximal for the Gaussian distribution with the same variance ${\overline{\sigma }}^2$, see [4]. Since the differential entropy of a Gaussian is $log\left(\overline{\sigma }\sqrt{2\pi e}\right)$, we deduce that the entropy of the GMM is upper bounded by

$$H\left(m\right)\le \frac{1}{2}log\left(2\pi e\right)+\frac{1}{2}log\sum _{i=1}^k{w}_i\left[{({\mu }_i-\overline{\mu })}^2+{\sigma }_i^2\right].$$

This upper bound can be easily generalized to arbitrary dimensionality. We get the following lemma:

Lemma 2.

The entropy of a d-variate GMM $m\left(x\right)={\sum }_{i=1}^k{w}_{i}p(x;{\mu }_i,{\mathsf{\Sigma }}_i)$ is upper bounded by $\frac{d}{2}log\left(2\pi e\right)+\frac{1}{2}logdet\mathsf{\Sigma }$, where $\mathsf{\Sigma }={\sum }_{i=1}^k{w}_i({\mu }_i{\mu }_i^{\top }+{\mathsf{\Sigma }}_i)-\left({\sum }_{i=1}^k{w}_i{\mu }_i\right)\left({\sum }_{i=1}^k{w}_i{\mu }_i^{\top }\right)$.

In general, exponential families have finite moments of any order [17]: In particular, we have $E\left[t\right(X\left)\right]=\nabla F\left(\theta \right)$ and $V\left[t\left(X\right)\right]={\nabla }^{2}F\left(\theta \right)$. For Gaussian distribution, we have the sufficient statistics $t\left(x\right)=(x,x^2)$ so that $E\left[t\right(X\left)\right]=\nabla F\left(\theta \right)$ yields the mean and variance from the log-normalizer. It is easy to generalize Lemma 2 to mixtures of exponential family distributions.

Note that this bound (called the Maximum Entropy Upper Bound in [13], MEUB) is tight when the GMM approximates a single Gaussian. It is fast to compute compared to the bound reported in [9] that uses Taylor’ s expansion of the log-sum of the mixture density.

A similar argument cannot be applied for a lower bound since a GMM with a given variance may have entropy tending to $-\infty$. For example, assume the 2-component mixture’s mean is zero, and that the variance approximates 1 by taking $m\left(x\right)=\frac{1}{2}G(x;-1,ϵ)+\frac{1}{2}G(x;1,ϵ)$ where G denotes the Gaussian density. Letting $ϵ\to 0$, we get the entropy tending to $-\infty$.

We remark that our log-sum-exp inequality technique yields a $log2$ additive approximation range in the case of a Gaussian mixture with two components. It thus generalizes the bounds reported in [7] to GMMs with arbitrary variances that are not necessarily equal.

To see the bound gap, we have

$$\begin{array}{c}-\sum _r{\int }_{I_r}m\left(x\right)\left(logk+log\underset{i}{max}{w}_i{p}_i\left(x\right)\right)\mathrm{d}x\le H\left(m\right)\hfill \\ \le -\sum _r{\int }_{I_r}m\left(x\right)max\left\{log\underset{i}{max}{w}_i{p}_i\left(x\right),logk+log\underset{i}{min}{w}_i{p}_i\left(x\right)\right\}\mathrm{d}x.\hfill \end{array}$$

(27)

Therefore the gap is at most

$$\begin{array}{cc}\hfill \Delta & =min\left\{\sum _r{\int }_{I_r}m\left(x\right)log\frac{{max}_i{w}_i{p}_i\left(x\right)}{{min}_i{w}_i{p}_i\left(x\right)}\mathrm{d}x,logk\right\}\hfill \\ & =min\left\{\sum _s\sum _r{\int }_{I_r}{w}_s{p}_s\left(x\right)log\frac{{max}_i{w}_i{p}_i\left(x\right)}{{min}_i{w}_i{p}_i\left(x\right)}\mathrm{d}x,logk\right\}.\hfill \end{array}$$

(28)

Thus to compute the gap error bound of the differential entropy, we need to integrate terms in the form

$$\int w_a{p}_a\left(x\right)log\frac{w_b{p}_b\left(x\right)}{w_c{p}_c\left(x\right)}\mathrm{d}x.$$

See Appendix B for a closed-form formula when dealing with exponential family components.

4. Bounding the α-Divergence

The α-divergence [15,32,33,34] between $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_i\left(x\right)$ and $m^{\prime }\left(x\right)={\sum }_{i=1}^{k^{\prime }}{w}_i^{\prime }{p}_i^{\prime }\left(x\right)$ is defined as

$$\begin{array}{c}\hfill D_{\alpha }\left(m:m^{\prime }\right)=\frac{1}{\alpha (1-\alpha )}\left(1-{\int }_{\mathcal{X}}m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }\mathrm{d}x\right),\end{array}$$

(29)

which clearly satisfies $D_{\alpha }\left(m:m^{\prime }\right)=D_{1-\alpha }\left(m^{\prime }:m\right)$. The α-divergence is a family of information divergences parametrized by $\alpha \in \mathbb{R}\backslash \{0,1\}$. Let $\alpha \to 1$, we get the KL divergence (see [35] for a proof):

$$\underset{\alpha \to 1}{lim}{D}_{\alpha }(m:m^{\prime })=\mathrm{KL}(m:m^{\prime })={\int }_{\mathcal{X}}m\left(x\right)log\frac{m\left(x\right)}{m^{\prime }\left(x\right)}\mathrm{d}x,$$

(30)

and $\alpha \to 0$ gives the reverse KL divergence:

$$\underset{\alpha \to 0}{lim}{D}_{\alpha }(m:m^{\prime })=\mathrm{KL}(m^{\prime }:m).$$

Other interesting values [33] include $\alpha =1/2$ (squared Hellinger distance), $\alpha =2$ (Pearson Chi-square distance), $\alpha =-1$ (Neyman Chi-square distance), etc. Notably, the Hellinger distance is a valid distance metric which satisfies non-negativity, symmetry, and the triangle inequality. In general, $D_{\alpha }(m:m^{\prime })$ only satisfies non-negativity so that $D_{\alpha }\left(m:m^{\prime }\right)\ge 0$ for any $m\left(x\right)$ and $m^{\prime }\left(x\right)$. It is neither symmetric nor admitting the triangle inequality. Minimization of α-divergences allows one to choose a trade-off between mode fitting and support fitting of the minimizer [36]. The minimizer of α-divergences including MLE as a special case has interesting connections with transcendental number theory [37].

To compute $D_{\alpha }\left(m:m^{\prime }\right)$ for given $m\left(x\right)$ and $m^{\prime }\left(x\right)$ reduces to evaluate the Hellinger integral [38,39]:

$$H_{\alpha }(m:m^{\prime })={\int }_{\mathcal{X}}m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }\mathrm{d}x,$$

(31)

which in general does not have a closed form, as it was known that the α-divergence of mixture models is not analytic [6]. Moreover, $H_{\alpha }(m:m^{\prime })$ may diverge making the α-divergence unbounded. Once $H_{\alpha }(m:m^{\prime })$ can be solved, the Rényi and Tsallis divergences [35] and in general Sharma–Mittal divergences [40] can be easily computed. Therefore the results presented here directly extend to those divergence families.

Similar to the case of KL divergence, the Monte Carlo stochastic estimation of $H_{\alpha }(m:m^{\prime })$ can be computed either as

$${\widehat{H}}_{\alpha }^n\left(m:m^{\prime }\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{m^{\prime }\left(x_i\right)}{m\left(x_i\right)}\right)}^{1-\alpha },$$

where $x_1,\dots ,x_n\sim m\left(x\right)$ are i.i.d. samples, or as

$${\widehat{H}}_{\alpha }^n\left(m:m^{\prime }\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{m\left(x_i\right)}{m^{\prime }\left(x_i\right)}\right)}^{\alpha },$$

where $x_1,\dots ,x_n\sim m^{\prime }\left(x\right)$ are i.i.d. In either case, it is consistent so that ${lim}_{n\to \infty }{\widehat{H}}_{\alpha }^n\left(m:m^{\prime }\right)=H_{\alpha }\left(m:m^{\prime }\right)$. However, MC estimation requires a large sample and does not guarantee deterministic bounds. The techniques described in [41] work in practice for very close distributions, and do not apply between mixture models. We will therefore derive combinatorial bounds for $H_{\alpha }(m:m^{\prime })$. The structure of this Section is parallel with Section 2 with necessary reformulations for a clear presentation.

4.1. Basic Bounds

For a pair of given $m\left(x\right)$ and $m^{\prime }\left(x\right)$, we only need to derive bounds of $H_{\alpha }(m:m^{\prime })$ in Equation (31) so that $L_{\alpha }(m:m^{\prime })\le H_{\alpha }(m:m^{\prime })\le U_{\alpha }(m:m^{\prime })$. Then the α-divergence $D_{\alpha }(m:m^{\prime })$ can be bounded by a linear transformation of the range $\left[L_{\alpha }(m:m^{\prime }),U_{\alpha }(m:m^{\prime })\right]$. In the following we always assume without loss of generality $\alpha \ge 1/2$. Otherwise we can bound $D_{\alpha }(m:m^{\prime })$ by considering equivalently the bounds of $D_{1-\alpha }(m^{\prime }:m)$.

Recall that in each elementary slab $I_r$, we have

$$\begin{array}{c}\hfill max\left\{k{w}_{ϵ\left(r\right)}{p}_{ϵ\left(r\right)}\left(x\right),w_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)\right\}\le m\left(x\right)\le k{w}_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right).\end{array}$$

(32)

Notice that $k{w}_{ϵ\left(r\right)}{p}_{ϵ\left(r\right)}\left(x\right)$, $w_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)$, and $k{w}_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)$ are all single component distributions up to a scaling coefficient. The general thinking is to bound the multi-component mixture $m\left(x\right)$ by single component distributions in each elementary interval, so that the integral in Equation (31) can be computed in a piecewise manner.

For the convenience of notation, we rewrite Equation (32) as

$$c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right)\le m\left(x\right)\le c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right),$$

(33)

where

$$c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right):=k{w}_{ϵ\left(r\right)}{p}_{ϵ\left(r\right)}\left(x\right)\mathrm{or}{w}_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right),$$

(34)

$$c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right):=k{w}_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right).$$

(35)

If $1/2\le \alpha <1$, then both $x^{\alpha }$ and $x^{1-\alpha }$ are monotonically increasing on ${\mathbb{R}}^{+}$. Therefore we have

$$\begin{array}{c}\hfill A_{\nu \left(r\right),{\nu }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right)\le {\int }_{I_r}m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }\mathrm{d}x\le A_{\delta \left(r\right),{\delta }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right),\end{array}$$

(36)

where

$$\begin{array}{c}\hfill A_{i,j}^{\alpha }\left(I\right)={\int }_I{\left(c_i{p}_i\left(x\right)\right)}^{\alpha }{\left(c_j^{\prime }{p}_j^{\prime }\left(x\right)\right)}^{1-\alpha }\mathrm{d}x,\end{array}$$

(37)

and I denotes an interval $I=(a,b)\subset \mathbb{R}$. The other case $\alpha >1$ is similar by noting that $x^{\alpha }$ and $x^{1-\alpha }$ are monotonically increasing and decreasing on ${\mathbb{R}}^{+}$, respectively. In conclusion, we obtain the following bounds of $H_{\alpha }(m:m^{\prime })$:

$$\begin{array}{cc}\hfill \mathrm{If}1/2\le \alpha <1,& L_{\alpha }(m:m^{\prime })=\sum _{r=1}^{\ell }{A}_{\nu \left(r\right),{\nu }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right),U_{\alpha }(m:m^{\prime })=\sum _{r=1}^{\ell }{A}_{\delta \left(r\right),{\delta }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right);\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill \mathrm{if}\alpha >1,& L_{\alpha }(m:m^{\prime })=\sum _{r=1}^{\ell }{A}_{\nu \left(r\right),{\delta }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right),U_{\alpha }(m:m^{\prime })=\sum _{r=1}^{\ell }{A}_{\delta \left(r\right),{\nu }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right).\hfill \end{array}$$

(39)

The remaining problem is to compute the definite integral $A_{i,j}^{\alpha }\left(I\right)$ in the above equations. Here we assume all mixture components are in the same exponential family so that $p_i\left(x\right)=p(x;{\theta }_i)=h\left(x\right)exp\left({\theta }_i^{\top }t\left(x\right)-F\left({\theta }_i\right)\right)$, where $h\left(x\right)$ is a base measure, $t\left(x\right)$ is a vector of sufficient statistics, and the function F is known as the cumulant generating function. Then it is straightforward from Equation (37) that

$$A_{i,j}^{\alpha }\left(I\right)=c_i^{\alpha }{\left(c_j^{\prime }\right)}^{1-\alpha }{\int }_{I}h\left(x\right)exp\left({\left(\alpha {\theta }_i+(1-\alpha ){\theta }_j^{\prime }\right)}^{\top }t\left(x\right)-\alpha F\left({\theta }_i\right)-(1-\alpha )F({\theta }_j^{\prime })\right)\mathrm{d}x.$$

(40)

If $1/2\le \alpha <1$, then $\overline{\theta }=\alpha {\theta }_i+(1-\alpha ){\theta }_j^{\prime }$ belongs to the natural parameter space ${\mathcal{M}}_{\theta }$. Therefore $A_{i,j}^{\alpha }\left(I\right)$ is bounded and can be computed from the CDF of $p(x;\overline{\theta })$ as

$$A_{i,j}^{\alpha }\left(I\right)=c_i^{\alpha }{\left(c_j^{\prime }\right)}^{1-\alpha }exp(F\left(\overline{\theta }\right)-\alpha F\left({\theta }_i\right)-(1-\alpha )F\left({\theta }_j^{\prime }\right)){\int }_{I}p\left(x;\overline{\theta }\right)\mathrm{d}x.$$

(41)

The other case $\alpha >1$ is more difficult: if $\overline{\theta }=\alpha {\theta }_i+(1-\alpha ){\theta }_j^{\prime }$ still lies in ${\mathcal{M}}_{\theta }$, then $A_{i,j}^{\alpha }\left(I\right)$ can be computed by Equation (41). Otherwise we try to solve it by a numerical integrator. This is not ideal as the integral may diverge, or our approximation may be too loose to conclude. We point the reader to [42] and Equations (61)–(69) in [35] for related analysis with more details. As computing $A_{i,j}^{\alpha }\left(I\right)$ only requires $O\left(1\right)$ time, the overall computational complexity (without considering the envelope computation) is $O\left(\ell \right)$.

4.2. Adaptive Bounds

This section derives the shape-dependent bounds which improve the basic bounds in Section 4.1. We can rewrite a mixture model $m\left(x\right)$ in a slab $I_r$ as

$$m\left(x\right)=w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)\left(1+\sum _{i\ne \zeta \left(r\right)}\frac{w_i{p}_i\left(x\right)}{w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)}\right),$$

(42)

where $w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)$ is a weighted component in $m\left(x\right)$ serving as a reference. We only discuss the case that the reference is chosen as the dominating component, i.e., $\zeta \left(r\right)=\delta \left(r\right)$. However it is worth to note that the proposed bounds do not depend on this particular choice. Therefore the ratio

$$\frac{w_i{p}_i\left(x\right)}{w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)}=\frac{w_i}{w_{\zeta \left(r\right)}}exp\left({\left({\theta }_i-{\theta }_{\zeta \left(r\right)}\right)}^{\top }t\left(x\right)-F\left({\theta }_i\right)+F\left({\theta }_{\zeta \left(r\right)}\right)\right)$$

(43)

can be bounded in a sub-range of $[0,1]$ by analyzing the extreme values of $t\left(x\right)$ in the slab $I_r$. This can be done because $t\left(x\right)$ usually consists of polynomial functions with finite critical points which can be solved easily. Correspondingly the function $\left(1+{\sum }_{i\ne \zeta \left(r\right)}\frac{w_i{p}_i\left(x\right)}{w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)}\right)$ in $I_r$ can be bounded in a subrange of $[1,k]$, denoted as $[{\omega }_{\zeta \left(r\right)}\left(I_r\right),{\mathsf{\Omega }}_{\zeta \left(r\right)}\left(I_r\right)]$. Hence

$$\begin{array}{c}\hfill {\omega }_{\zeta \left(r\right)}\left(I_r\right)w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right)\le m\left(x\right)\le {\mathsf{\Omega }}_{\zeta \left(r\right)}\left(I_r\right)w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right).\end{array}$$

(44)

This forms better bounds of $m\left(x\right)$ than Equation (32) because each component in the slab $I_r$ is analyzed more accurately. Therefore, we refine the fundamental bounds of $m\left(x\right)$ by replacing the Equations (34) and (35) with

$$\begin{array}{c}\hfill c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right):={\omega }_{\zeta \left(r\right)}\left(I_r\right)w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right),\end{array}$$

(45)

$$\begin{array}{c}\hfill c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right):={\mathsf{\Omega }}_{\zeta \left(r\right)}\left(I_r\right)w_{\zeta \left(r\right)}{p}_{\zeta \left(r\right)}\left(x\right).\end{array}$$

(46)

Then, the improved bounds of $H_{\alpha }$ are given by Equations (38) and (39) according to the above replaced definition of $c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right)$ and $c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)$.

To evaluate ${\omega }_{\zeta \left(r\right)}\left(I_r\right)$ and ${\mathsf{\Omega }}_{\zeta \left(r\right)}\left(I_r\right)$ requires iterating through all components in each slab. Therefore the computational complexity is increased to $O\left(\ell (k+k^{\prime })\right)$.

4.3. Variance-Reduced Bounds

This section further improves the proposed bounds based on variance reduction [43]. By assumption, $\alpha \ge 1/2$, then $m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }$ is more similar to $m\left(x\right)$ rather than $m^{\prime }\left(x\right)$. The ratio $m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }/m\left(x\right)$ is likely to have a small variance when x varies inside a slab $I_r$, especially when α is close to 1. We will therefore bound this ratio term in

$$\begin{array}{c}{\int }_{I_r}m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }\mathrm{d}x={\int }_{I_r}m\left(x\right)\left(\frac{m{\left(x\right)}^{\alpha }{m}^{\prime }{\left(x\right)}^{1-\alpha }}{m\left(x\right)}\right)\mathrm{d}x=\sum _{i=1}^k{\int }_{I_r}{w}_i{p}_i\left(x\right){\left(\frac{m^{\prime }\left(x\right)}{m\left(x\right)}\right)}^{1-\alpha }\mathrm{d}x.\hfill \end{array}$$

(47)

No matter $\alpha <1$ or $\alpha >1$, the function $x^{1-\alpha }$ must be monotonic on ${\mathbb{R}}^{+}$. In each slab $I_r$, ${\left(m^{\prime }\left(x\right)/m\left(x\right)\right)}^{1-\alpha }$ ranges between these two functions:

$$\begin{array}{c}\hfill {\left(\frac{c_{{\nu }^{\prime }\left(r\right)}^{\prime }{p}_{{\nu }^{\prime }\left(r\right)}^{\prime }\left(x\right)}{c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)}\right)}^{1-\alpha }\mathrm{and}{\left(\frac{c_{{\delta }^{\prime }\left(r\right)}^{\prime }{p}_{{\delta }^{\prime }\left(r\right)}^{\prime }\left(x\right)}{c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right)}\right)}^{1-\alpha },\end{array}$$

(48)

where $c_{\nu \left(r\right)}{p}_{\nu \left(r\right)}\left(x\right)$, $c_{\delta \left(r\right)}{p}_{\delta \left(r\right)}\left(x\right)$, $c_{{\nu }^{\prime }\left(r\right)}^{\prime }{p}_{{\nu }^{\prime }\left(r\right)}^{\prime }\left(x\right)$ and $c_{{\delta }^{\prime }\left(r\right)}^{\prime }{p}_{{\delta }^{\prime }\left(r\right)}^{\prime }\left(x\right)$ are defined in Equations (45) and (46). Similar to the definition of $A_{i,j}^{\alpha }\left(I\right)$ in Equation (37), we define

$$\begin{array}{c}\hfill B_{i,j,l}^{\alpha }\left(I\right)={\int }_I{w}_i{p}_i\left(x\right){\left(\frac{c_l^{\prime }{p}_l^{\prime }\left(x\right)}{c_j{p}_j\left(x\right)}\right)}^{1-\alpha }\mathrm{d}x.\end{array}$$

(49)

Therefore we have,

$$\begin{array}{c}{L}_{\alpha }(m:m^{\prime })=min\mathcal{S},U_{\alpha }(m:m^{\prime })=max\mathcal{S},\hfill \\ \mathcal{S}=\left\{\sum _{r=1}^{\ell }\sum _{i=1}^k{B}_{i,\delta \left(r\right),{\nu }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right),\sum _{r=1}^{\ell }\sum _{i=1}^k{B}_{i,\nu \left(r\right),{\delta }^{\prime }\left(r\right)}^{\alpha }\left(I_r\right)\right\}.\hfill \end{array}$$

(50)

The remaining problem is to evaluate $B_{i,j,l}^{\alpha }\left(I\right)$ in Equation (49). Similar to Section 4.1, assuming the components are in the same exponential family with respect to the natural parameters θ, we get

$$\begin{array}{cc}\hfill B_{i,j,l}^{\alpha }\left(I\right)=& w_i\frac{c_l^{\prime 1-\alpha }}{c_j^{1-\alpha }}exp\left(F\left(\overline{\theta }\right)-F\left({\theta }_i\right)-(1-\alpha )F\left({\theta }_l^{\prime }\right)+(1-\alpha )F\left({\theta }_j\right)\right){\int }_{I}p(x;\overline{\theta })\mathrm{d}x.\hfill \end{array}$$

(51)

If $\overline{\theta }={\theta }_i+(1-\alpha ){\theta }_l^{\prime }-(1-\alpha ){\theta }_j$ is in the natural parameter space, $B_{i,j,l}^{\alpha }\left(I\right)$ can be computed from the CDF of $p(x;\overline{\theta })$; otherwise $B_{i,j,l}^{\alpha }\left(I\right)$ can be numerically integrated by its definition in Equation (49). The computational complexity is the same as the bounds in Section 4.2, i.e., $O\left(\ell \left(k+k^{\prime }\right)\right)$.

We have introduced three pairs of deterministic lower and upper bounds that enclose the true value of α-divergence between univariate mixture models. Thus the gap between the upper and lower bounds provides the additive approximation factor of the bounds. We conclude by emphasizing that the presented methodology can be easily generalized to other divergences [35,40] relying on Hellinger-type integrals $H_{\alpha ,\beta }(p:q)=\int p{\left(x\right)}^{\alpha }q{\left(x\right)}^{\beta }\mathrm{d}x$ like the γ-divergence [44] as well as entropy measures [45].

5. Lower Bounds of the $\mathit{f}$-Divergence

The f-divergence between two distributions $m\left(x\right)$ and $m^{\prime }\left(x\right)$ (not necessarily mixtures) is defined for a convex generator f by:

$$D_f(m:m^{\prime })=\int m\left(x\right)f\left(\frac{m^{\prime }\left(x\right)}{m\left(x\right)}\right)\mathrm{d}x.$$

If $f\left(x\right)=-logx$, then $D_f(m:m^{\prime })=\mathrm{KL}(m:m^{\prime })$.

Let us partition the support $\mathcal{X}={\uplus }_{r=1}^{\ell }{I}_r$ arbitrarily into elementary ranges, which do not necessarily correspond to the envelopes. Denote by $M_I$ the probability mass of a mixture $m\left(x\right)$ in the range I: $M_I={\int }_{I}m\left(x\right)dx$. Then

$$D_f(m:m^{\prime })=\sum _{r=1}^{\ell }{M}_{I_r}{\int }_{I_r}\frac{m\left(x\right)}{M_{I_r}}f\left(\frac{m^{\prime }\left(x\right)}{m\left(x\right)}\right)\mathrm{d}x.$$

Note that in range $I_r$, $\frac{m\left(x\right)}{M_{I_r}}$ is a unit weight distribution. Thus by Jensen’s inequality $f\left(E\right[X\left]\right)\le E\left[f\right(X\left)\right]$, we get

$$D_f(m:m^{\prime })\ge \sum _{r=1}^{\ell }{M}_{I_r}f\left({\int }_{I_r}\frac{m\left(x\right)}{M_{I_r}}\frac{m^{\prime }\left(x\right)}{m\left(x\right)}\mathrm{d}x\right)=\sum _{r=1}^{\ell }{M}_{I_r}f\left(\frac{M_{I_r}^{\prime }}{M_{I_r}}\right).$$

(52)

Notice that the RHS of Equation (52) is the f-divergence between $(M_{I_1},\cdots ,M_{I_{\ell }})$ and $(M_{I_1}^{\prime },\cdots ,M_{I_{\ell }}^{\prime })$, denoted by $D_f^{\mathcal{I}}(m:m^{\prime })$. In the special case that $\ell =1$ and $I_1=\mathcal{X}$, the above Equation (52) turns out to be the usual Gibbs’ inequality: $D_f(m:m^{\prime })\ge f\left(1\right)$, and Csiszár generator is chosen so that $f\left(1\right)=0$. In conclusion, for a fixed (coarse-grained) countable partition of $\mathcal{X}$, we recover the well-know information monotonicity [46] of the f-divergences:

$$D_f(m:m^{\prime })\ge D_f^{\mathcal{I}}(m:m^{\prime })\ge 0.$$

In practice, we get closed-form lower bounds when $M_I={\int }_a^{b}m\left(x\right)\mathrm{d}x=\mathrm{\Phi }\left(b\right)-\mathrm{\Phi }\left(a\right)$ is available in closed-form, where $\mathrm{\Phi }(·)$ denote the CDF. In particular, if $m\left(x\right)$ is a mixture model, then its CDF can be computed by linearly combining the CDFs of its components.

To wrap up, we have proved that coarse-graining by making a finite partition of the support $\mathcal{X}$ yields a lower bound on the f-divergence by virtue of the information monotonicity. Therefore, instead of doing Monte Carlo stochastic integration:

$${\widehat{D}}_f^n(m:m^{\prime })=\frac{1}{n}\sum _{i=1}^{n}f\left(\frac{m^{\prime }\left(x_i\right)}{m\left(x_i\right)}\right),$$

with $x_1,\dots ,x_n{\sim }_{\text{i.i.d.}}m\left(x\right)$, it could be better to sort those n samples and consider the coarse-grained partition:

$$\mathcal{I}=(-\infty ,x_{\left(1\right)}]\cup \left({\uplus }_{i=1}^{n-1}(x_{\left(i\right)},x_{(i+1)}]\right)\cup (x_{\left(n\right)},+\infty )$$

to get a guaranteed lower bound on the f-divergence. We will call this bound CGQLB for Coarse Graining Quantization Lower Bound.

Given a budget of n splitting points on the range $\mathcal{X}$, it would be interesting to find the best n points that maximize the lower bound $D_f^{\mathcal{I}}(m:m^{\prime })$. This is ongoing research.

6. Experiments

We perform an empirical study to verify our theoretical bounds. We simulate four pairs of mixture models $\{({\mathtt{EMM}}_1,{\mathtt{EMM}}_2),({\mathtt{RMM}}_1,{\mathtt{RMM}}_2),({\mathtt{GMM}}_1,{\mathtt{GMM}}_2),({\mathtt{GaMM}}_1,{\mathtt{GaMM}}_2)\}$ as the test subjects. The component type is implied by the model name, where $\mathtt{GaMM}$ stands for Gamma mixtures. The components of each mixture model are given as follows.

• ${\mathtt{EMM}}_1$’s components, in the form $({\lambda }_i,w_i)$, are given by $(0.1,1/3)$, $(0.5,1/3)$, $(1,1/3)$; ${\mathtt{EMM}}_2$’s components are $(2,0.2)$, $(10,0.4)$, $(20,0.4)$.
• ${\mathtt{RMM}}_1$’s components, in the form $({\sigma }_i,w_i)$, are given by $(0.5,1/3)$, $(2,1/3)$, $(10,1/3)$; ${\mathtt{RMM}}_2$ consists of $(5,0.25)$, $(60,0.25)$, $(100,0.5)$.
• ${\mathtt{GMM}}_1$’s components, in the form $({\mu }_i,{\sigma }_i,w_i)$, are $(-5,1,0.05)$, $(-2,0.5,0.1)$, $(5,0.3,0.2)$, $(10,0.5,0.2)$, $(15,0.4,0.05)$, $(25,0.5,0.3)$, $(30,2,0.1)$; ${\mathtt{GMM}}_2$ consists of $(-16,0.5,0.1)$, $(-12,0.2,0.1)$, $(-8,0.5,0.1)$, $(-4,0.2,0.1)$, $(0,0.5,0.2)$, $(4,0.2,0.1)$, $(8,0.5,0.1)$, $(12,0.2,0.1)$, $(16,0.5,0.1)$.
• ${\mathtt{GaMM}}_1$’s components, in the form $(k_i,{\lambda }_i,w_i)$, are $(2,0.5,1/3)$, $(2,2,1/3)$, $(2,4,1/3)$; ${\mathtt{GaMM}}_2$ consists of $(2,5,1/3)$, $(2,8,1/3)$, $(2,10,1/3)$.

We compare the proposed bounds with Monte Carlo estimation with different sample sizes in the range $\{{10}^2,{10}^3,{10}^4,{10}^5\}$. For each sample size configuration, we report the 0.95 confidence interval by Monte Carlo estimation using the corresponding number of samples. Figure 2a–d shows the input signals as well as the estimation results, where the proposed bounds CELB, CEUB, CEALB, CEAUB, CGQLB are presented as horizontal lines, and the Monte Carlo estimations over different sample sizes are presented as error bars. We can loosely consider the average Monte Carlo output with the largest sample size (${10}^5$) as the underlying truth, which is clearly inside our bounds. This serves as an empirical justification on the correctness of the bounds.

A key observation is that the bounds can be very tight, especially when the underlying KL divergence has a large magnitude, e.g., $\mathrm{KL}({\mathtt{RMM}}_2:{\mathtt{RMM}}_1)$. This is because the gap between the lower and upper bounds is always guaranteed to be within $logk+log{k}^{\prime }$. Because KL is unbounded [4], in the general case two mixture models may have a large KL. Then our approximation gap is relatively very small. On the other hand, we also observed that the bounds in certain cases, e.g., $\mathrm{KL}({\mathtt{EMM}}_2:{\mathtt{EMM}}_1)$, are not as tight as the other cases. When the underlying KL is small, the bounds are not as informative as the general case.

Comparatively, there is a significant improvement of the shape-dependent bounds (CEALB and CEAUB) over the combinatorial bounds (CELB and CEUB). In all investigated cases, the adaptive bounds can roughly shrink the gap by half of its original size at the cost of additional computation.

Note that, the bounds are accurate and must contain the true value. Monte Carlo estimation gives no guarantee on where the true value is. For example, in estimating $\mathrm{KL}({\mathtt{GMM}}_1:{\mathtt{GMM}}_2)$, Monte Carlo estimation based on ${10}^4$ samples can go beyond our bounds! It therefore suffers from a larger estimation error.

CGQLB as a simple-to-implement technique shows surprising good performance in several cases, e.g., $\mathrm{KL}({\mathtt{RMM}}_1,{\mathtt{RMM}}_2)$. Although it requires a large number of samples, we can observe that increasing sample size has limited effect on improving this bound. Therefore, in practice, one may intersect the range defined by CEALB and CEAUB with the range defined by CGQLB with a small sample size (e.g., 100) to get better bounds.

We simulates a set of Gaussian mixture models besides the above ${\mathtt{GMM}}_1$ and ${\mathtt{GMM}}_2$. Figure 3 shows the GMM densities as well as their differential entropy. A detailed explanation of the components of each GMM model is omitted for brevity.

The key observation is that CEUB (CEAUB) is very tight in most of the investigated cases. This is because that the upper envelope that is used to compute CEUB (CEAUB) gives a very good estimation of the input signal.

Notice that MEUB only gives an upper bound of the differential entropy as discussed in Section 3. In general the proposed bounds are tighter than MEUB. However, this is not the case when the mixture components are merged together and approximate one single Gaussian (and therefore its entropy can be well approximated by the Gaussian entropy), as shown in the last line of Figure 3.

For α-divergence, the bounds introduced in Section 4.1, Section 4.2 and Section 4.3 are denoted as “Basic”, “Adaptive” and “VR”, respectively. Figure 4 visualizes these GMMs and plots the estimations of their α-divergences against α. The red lines mean the upper envelope. The dashed vertical lines mean the elementary intervals. The components of ${\mathtt{GMM}}_{\mathtt{1}}$ and ${\mathtt{GMM}}_{\mathtt{2}}$ are more separated than ${\mathtt{GMM}}_{\mathtt{3}}$ and ${\mathtt{GMM}}_{\mathtt{4}}$. Therefore these two pairs present different cases. For a clear presentation, only VR (which is expected to be better than Basic and Adaptive) is shown. We can see that, visually in the big scale, VR tightly surrounds the true value.

For a more quantitative comparison,

Table 1. The estimated $D_{\alpha }$ and its bounds. The 95% confidence interval is shown for MC.

	α	MC(${10}^2$)	MC(${10}^3$)	MC(${10}^4$)	Basic		Adaptive		VR
					L	U	L	U	L	U
${\mathtt{GMM}}_{\mathtt{1}}\&{\mathtt{GMM}}_{\mathtt{2}}$	0	$15.96\pm 3.9$	$12.30\pm 1.0$	$13.63\pm 0.3$	11.75	15.89	12.96	14.63
	$0.01$	$13.36\pm 2.9$	$10.63\pm 0.8$	$11.66\pm 0.3$	−700.50	11.73	−77.33	11.73	11.40	12.27
	$0.5$	$3.57\pm 0.3$	$3.47\pm 0.1$	$3.47\pm 0.07$	−0.60	3.42	3.01	3.42	3.17	3.51
	$0.99$	$40.04\pm 7.7$	$37.22\pm 2.3$	$38.58\pm 0.8$	−333.90	39.04	5.36	38.98	38.28	38.96
	1	$104.01\pm 28$	$84.96\pm 7.2$	$92.57\pm 2.5$	91.44	95.59	92.76	94.41
${\mathtt{GMM}}_{\mathtt{3}}\&{\mathtt{GMM}}_{\mathtt{4}}$	0	$0.71\pm 0.2$	$0.63\pm 0.07$	$0.62\pm 0.02$	0.00	1.76	0.00	1.16
	$0.01$	$0.71\pm 0.2$	$0.63\pm 0.07$	$0.62\pm 0.02$	−179.13	7.63	−38.74	4.96	0.29	1.00
	$0.5$	$0.82\pm 0.3$	$0.57\pm 0.1$	$0.62\pm 0.04$	−5.23	0.93	−0.71	0.85	−0.18	1.19
	$0.99$	$0.79\pm 0.3$	$0.76\pm 0.1$	$0.80\pm 0.03$	−165.72	12.10	−59.76	9.11	0.37	1.28
	1	$0.80\pm 0.3$	$0.77\pm 0.1$	$0.81\pm 0.03$	0.00	1.82	0.31	1.40

shows the estimated α-divergence by MC, Basic, Adaptive, and VR. As $D_{\alpha }$ is defined on $\mathbb{R}\backslash \{0,1\}$, the KL bounds CE(A)LB and CE(A)UB are presented for $\alpha =0$ or 1. Overall, we have the following order of gap size: Basic > Adaptive > VR, and VR is recommended in general for bounding α-divergences. There are certain cases that the upper VR bound is looser than Adaptive. In practice one can compute the intersection of these bounds as well as the trivial bound $D_{\alpha }(m:m^{\prime })\ge 0$ to get the best estimation.

Note the similarity between KL in Equation (30) and the expression in Equation (47). We give without a formal analysis that: CEAL(U)B is equivalent to VR at the limit $\alpha \to 0$ or $\alpha \to 1$. Experimentally as we slowly set $\alpha \to 1$, we can see that VR is consistent with CEAL(U)B.

7. Concluding Remarks and Perspectives

We have presented a fast versatile method to compute bounds on the Kullback–Leibler divergence between mixtures by building algorithmic formulae. We reported on our experiments for various mixture models in the exponential family. For univariate GMMs, we get a guaranteed bound of the KL divergence of two mixtures m and $m^{\prime }$ with k and $k^{\prime }$ components within an additive approximation factor of $logk+log{k}^{\prime }$ in $O\left((k+k^{\prime })log(k+k^{\prime })\right)$-time. Therefore, the larger the KL divergence, the better the bound when considering a multiplicative $(1+\alpha )$-approximation factor, since $\alpha =\frac{logk+log{k}^{\prime }}{\mathrm{KL}(m:m^{\prime })}$. The adaptive bounds are guaranteed to yield better bounds at the expense of computing potentially $O\left(k^2+{\left(k^{\prime }\right)}^2\right)$ intersection points of pairwise weighted components.

Our technique also yields the bound for the Jeffreys divergence (the symmetrized KL divergence: $J(m,m^{\prime })=\mathrm{KL}(m:m^{\prime })+\mathrm{KL}(m^{\prime }:m)$) and the Jensen–Shannon divergence [47] (JS):

$$\mathrm{JS}(m,m^{\prime })=\frac{1}{2}\left(\mathrm{KL}\left(m:\frac{m+m^{\prime }}{2}\right)+\mathrm{KL}\left(m^{\prime }:\frac{m+m^{\prime }}{2}\right)\right),$$

since $\frac{m+m^{\prime }}{2}$ is a mixture model with $k+k^{\prime }$ components. One advantage of this statistical distance is that it is symmetric, always bounded by $log2$, and its square root yields a metric distance [48]. The log-sum-exp inequalities may also be used to compute some Rényi divergences [35]:

$$R_{\alpha }(m,p)=\frac{1}{\alpha -1}log\left(\int m{\left(x\right)}^{\alpha }p{\left(x\right)}^{1-\alpha }\right)\mathrm{d}x,$$

when α is an integer, $m\left(x\right)$ a mixture and $p\left(x\right)$ a single (component) distribution. Getting fast guaranteed tight bounds on statistical distances between mixtures opens many avenues. For example, we may consider building hierarchical mixture models by merging iteratively two mixture components so that those pairs of components are chosen so that the KL distance between the full mixture and the simplified mixture is minimized.

In order to be useful, our technique is unfortunately limited to univariate mixtures: indeed, in higher dimensions, we can still compute the maximization diagram of weighted components (an additively weighted Bregman–Voronoi diagram [49,50] for components belonging to the same exponential family). However, it becomes more complex to compute in the elementary Voronoi cells V, the functions $C_{i,j}\left(V\right)$ and $M_i\left(V\right)$ (in 1D, the Voronoi cells are segments). We may obtain hybrid algorithms by approximating or estimating these functions. In 2D, it is thus possible to obtain lower and upper bounds on the Mutual Information [51] (MI) when the joint distribution $m(x,y)$ is a 2D mixture of Gaussians:

$$I(M;M^{\prime })=\int m(x,y)log\frac{m(x,y)}{m\left(x\right)m^{\prime }\left(y\right)}\mathrm{d}x\mathrm{d}y.$$

Indeed, the marginal distributions $m\left(x\right)$ and $m^{\prime }\left(y\right)$ are univariate Gaussian mixtures.

A Python code implementing those computational-geometric methods for reproducible research is available online [52].