Fast Approximations of the Jeffreys Divergence between Univariate Gaussian Mixtures via Mixture Conversions to Exponential-Polynomial Distributions

Abstract

The Jeffreys divergence is a renown arithmetic symmetrization of the oriented Kullback–Leibler divergence broadly used in information sciences. Since the Jeffreys divergence between Gaussian mixture models is not available in closed-form, various techniques with advantages and disadvantages have been proposed in the literature to either estimate, approximate, or lower and upper bound this divergence. In this paper, we propose a simple yet fast heuristic to approximate the Jeffreys divergence between two univariate Gaussian mixtures with arbitrary number of components. Our heuristic relies on converting the mixtures into pairs of dually parameterized probability densities belonging to an exponential-polynomial family. To measure with a closed-form formula the goodness of fit between a Gaussian mixture and an exponential-polynomial density approximating it, we generalize the Hyvärinen divergence to $\alpha$-Hyvärinen divergences. In particular, the 2-Hyvärinen divergence allows us to perform model selection by choosing the order of the exponential-polynomial densities used to approximate the mixtures. We experimentally demonstrate that our heuristic to approximate the Jeffreys divergence between mixtures improves over the computational time of stochastic Monte Carlo estimations by several orders of magnitude while approximating the Jeffreys divergence reasonably well, especially when the mixtures have a very small number of modes.

1. Introduction

1.1. Statistical Mixtures and Statistical Divergences

We consider the problem of approximating the Jeffreys divergence [1] between two finite univariate continuous mixture models [2] $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)$ with continuous component distributions $p_i$’s and $p_i^{″}$s defined on a coinciding support $\mathcal{X}\subset \mathbb{R}$. The mixtures $m\left(x\right)$ and $m^{\prime }\left(x\right)$ may have a different number of components (i.e., $k\ne k^{\prime }$). Historically, Pearson [3] first considered a univariate Gaussian mixture of two components for modeling the distribution of the ratio of forehead breadth to body length of a thousand crabs in 1894 (Pearson obtained a unimodal mixture).

Although our work applies to any continuous mixtures of an exponential family (e.g., Rayleigh mixtures [4] with restricted support $\mathcal{X}={\mathbb{R}}_{+}$), we explain our method for the most prominent family of mixtures encountered in practice: the Gaussian mixture models or GMMs for short. In the remainder, a univariate GMM $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_{{\mu }_i,{\sigma }_i}\left(x\right)$ with k Gaussian components

$$p_i\left(x\right)=p_{{\mu }_i,{\sigma }_i}\left(x\right):=\frac{1}{{\sigma }_i\sqrt{2\pi }}exp\left(-\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}\right),$$

is called a k-GMM.

The Kullback–Leibler divergence (KLD) [5,6] $D_{\mathrm{KL}}[m:m^{\prime }]$ between two mixtures m and $m^{\prime }$ is:

$$D_{\mathrm{KL}}[m:m^{\prime }]:={\int }_{\mathcal{X}}m\left(x\right)log\left(\frac{m\left(x\right)}{m^{\prime }\left(x\right)}\right)\mathrm{d}x.$$

(1)

The KLD is an oriented divergence since $D_{\mathrm{KL}}[m:m^{\prime }]\ne D_{\mathrm{KL}}[m^{\prime }:m]$.

The Jeffreys divergence (JD) [1] $D_J[m,m^{\prime }]$ is the arithmetic symmetrization of the forward KLD and the reverse KLDs:

$$\begin{array}{ccc}\hfill D_J[m,m^{\prime }]& :=& D_{\mathrm{KL}}[m:m^{\prime }]+D_{\mathrm{KL}}[m^{\prime }:m],\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& =& {\int }_{\mathcal{X}}(m\left(x\right)-m^{\prime }\left(x\right))log\left(\frac{m\left(x\right)}{m^{\prime }\left(x\right)}\right)\mathrm{d}x.\hfill \end{array}$$

(3)

The JD is a symmetric divergence: $D_J[m,m^{\prime }]=D_J[m^{\prime },m]$. In the literature, the Jeffreys divergence [7] has also been called the J-divergence [8,9], the symmetric Kullback–Leibler divergence [10] and sometimes the symmetrical Kullback–Leibler divergence [11,12]. In general, it is provably hard to calculate the definite integral of the KLD between two continuous mixtures in closed-form: For example, the KLD between two GMMs has been shown to be non-analytic [13]. Thus, in practice, when calculating the JD between two GMMs, one can either approximate [14,15], estimate [16], or bound [17,18] the KLD between mixtures. Another approach to bypass the computational intractability of calculating the KLD between mixtures consists of designing new types of divergences that admit closed-form expressions for mixtures. See, for example, the Cauchy–Schwarz divergence [19] or the total square divergence [20] (a total Bregman divergence) that admit the closed-form formula when handling GMMs. The total square divergence [20] is invariant to rigid transformations and provably robust to outliers in clustering applications.

In practice, to estimate the KLD between mixtures, one uses the following Monte Carlo (MC) estimator:

$${\widehat{D}}_{\mathrm{KL}}^{{\mathcal{S}}_s}[m:m^{\prime }]:=\frac{1}{s}\sum _{i=1}^s\left(log\left(\frac{m\left(x_i\right)}{m^{\prime }\left(x_i\right)}\right)+\frac{m^{\prime }\left(x_i\right)}{m\left(x_i\right)}-1\right)\ge 0,$$

where ${\mathcal{S}}_s=\{x_1,\dots ,x_s\}$ is s independent and identically distributed (i.i.d.) samples from $m\left(x\right)$. This MC estimator is by construction always non-negative and therefore consistent. That is, we have ${lim}_{s\to \infty }{\widehat{D}}_{\mathrm{KL}}^{{\mathcal{S}}_s}[m:m^{\prime }]=D_{\mathrm{KL}}[m:m^{\prime }]$ under mild conditions [21].

Similarly, we estimate the Jeffreys divergence via MC sampling as follows:

$${\widehat{D}}_J^{{\mathcal{S}}_s}[m,m^{\prime }]:=\frac{1}{s}\sum _{i=1}^{s}2\frac{(m\left(x_i\right)-m^{\prime }\left(x_i\right))}{m\left(x_i\right)+m^{\prime }\left(x_i\right)}log\left(\frac{m\left(x_i\right)}{m^{\prime }\left(x_i\right)}\right)\ge 0,$$

(4)

where ${\mathcal{S}}_s=\{x_1,\dots ,x_s\}$ are s i.i.d. samples from the “middle mixture” $m_{12}\left(x\right):=\frac{1}{2}(m\left(x\right)+m^{\prime }\left(x\right))$. By choosing the middle mixture $m_{12}\left(x\right)$ for sampling, we ensure that we keep the symmetric property of the JD (i.e., ${\widehat{D}}_J^{{\mathcal{S}}_s}[m,m^{\prime }]={\widehat{D}}_J^{{\mathcal{S}}_s}[m^{\prime },m]$), and we also have consistency under mild conditions [21]: ${lim}_{s\to \infty }{\widehat{D}}_J^{{\mathcal{S}}_s}[m,m^{\prime }]=D_J[m,m^{\prime }]$. The time complexity to stochastically estimate the JD is $\tilde{O}\left((k+k^{\prime })s\right)$, with s typically ranging from ${10}^4$ to ${10}^6$ in applications. Notice that the number of components of a mixture can be very large (e.g., $k=O\left(n\right)$ for n input data when using Kernel Density Estimators [2]). KDEs may also have a large number of components and may potentially exhibit many spurious modes visualized as small bumps when plotting the densities.

1.2. Jeffreys Divergence between Densities of an Exponential Family

We consider approximating the JD by converting continuous mixtures into densities of exponential families [22]. A continuous exponential family (EF) ${\mathcal{E}}_t$ of order D is defined as a family of probability density functions with support $\mathcal{X}$ and the probability density function:

$${\mathcal{E}}_t:=\left\{p_{\theta }\left(x\right):=exp\left(\sum _{i=1}^D{\theta }_i{t}_i\left(x\right)-F\left(\theta \right)\right):\theta \in \Theta \right\},$$

where $F\left(\theta \right)$ is called the log-normalizer, which ensures the normalization of $p_{\theta }\left(x\right)$ (i.e., ${\int }_{\mathcal{X}}{p}_{\theta }\left(x\right)\mathrm{d}x=1$):

$$F\left(\theta \right)=log\left({\int }_{\mathcal{X}}exp\left(\sum _{i=1}^D{\theta }_i{t}_i\left(x\right)\right)\mathrm{d}x\right).$$

Parameter $\theta \in \Theta \subset {\mathbb{R}}^D$ is called the natural parameter, and the functions $t_1\left(x\right)$, …, $t_D\left(x\right)$ are called the sufficient statistics [22]. Let $\Theta$ denote the natural parameter space: $\Theta :=\{\theta :F(\theta )<\infty \}$, an open convex domain for regular exponential families [22]. The exponential family is said to be minimal when the functions $1,t_1\left(x\right),\dots ,t_D\left(x\right)$ are linearly independent.

It is well-known that one can bypass the definite integral calculation of the KLD when the probability density functions $p_{\theta }$ and $p_{{\theta }^{\prime }}$ belong to the same exponential family [23,24]:

$$D_{\mathrm{KL}}[p_{\theta }:p_{\theta }]=B_F({\theta }^{\prime }:\theta ),$$

where $B_F({\theta }_2:{\theta }_1)$ is the Bregman divergence induced by the log-normalizer, a strictly convex real-analytic function [22]. The Bregman divergence [25] between two parameters ${\theta }_1$ and ${\theta }_2$ for a strictly convex and smooth generator F is defined by:

$$B_F({\theta }_1:{\theta }_2):=F\left({\theta }_1\right)-F\left({\theta }_2\right)-{({\theta }_1-{\theta }_2)}^{\top }\nabla F\left({\theta }_2\right).$$

(5)

Thus, the Jeffreys divergence between two pdfs $p_{\theta }$ and $p_{{\theta }^{\prime }}$ belonging to the same exponential family is a symmetrized Bregman divergence [26]:

$$\begin{array}{ccc}\hfill D_J[p_{\theta }:p_{\theta }]& =& B_F({\theta }^{\prime }:\theta )+B_F(\theta :{\theta }^{\prime }),\hfill \\ & =& {({\theta }^{\prime }-\theta )}^{\top }(\nabla F\left({\theta }^{\prime }\right)-\nabla F\left(\theta \right)).\hfill \end{array}$$

Let $F^{*}\left(\eta \right)$ denote the Legendre–Fenchel convex conjugate of $F\left(\theta \right)$:

$$F^{*}\left(\eta \right):=\underset{\theta \in \Theta }{sup}\{{\theta }^{\top }\eta -F\left(\theta \right)\}.$$

(6)

The Legendre transform ensures that $\eta =\nabla F\left(\theta \right)$ and $\theta =\nabla F^{*}\left(\eta \right)$, and the Jeffreys divergence between two pdfs $p_{\theta }$ and $p_{{\theta }^{\prime }}$ belonging to the same exponential family is:

$$D_J[p_{\theta }:p_{\theta }]={({\theta }^{\prime }-\theta )}^{\top }({\eta }^{\prime }-\eta ).$$

(7)

Notice that the log-normalizer $F\left(\theta \right)$ does not appear explicitly in the above formula.

1.3. A Simple Approximation Heuristic

Densities $p_{\theta }$ of an exponential family admit a dual parameterization [22]: $\eta =\eta \left(\theta \right):=E_{p_{\theta }}\left[t\left(x\right)\right]=\nabla F\left(\theta \right)$, called the moment parameterization (or mean parameterization). Let H denote the moment parameter space. Let us use the subscript and superscript notations to emphasize the coordinate system used to index a density: In our notation, we thus write $p_{\theta }\left(x\right)=p^{\eta }\left(x\right)$.

In view of Equation (7), our method to approximate the Jeffreys divergence between mixtures m and $m^{\prime }$ consists of first converting those mixtures m and $m^{\prime }$ into pairs of polynomial exponential densities (PEDs) in Section 2. To convert a mixture $m\left(x\right)$ into a pair $(p_{{\overline{\theta }}_1},p^{{\overline{\eta }}_2})$ dually parameterized (but not dual because ${\overline{\eta }}_2\ne \nabla F\left({\overline{\theta }}_1\right)$), we shall consider “integral extensions” (or information projections) of the Maximum Likelihood Estimator [22] (MLE estimates in the moment parameter space $H=\{\nabla F(\theta ):\theta \in \Theta \}$) and of the Score Matching Estimator [27] (SME estimates in the natural parameter space $\Theta =\{\nabla F^{*}\left(\eta \right):\eta \in H\}$).

We shall consider polynomial exponential families [28] (PEFs) also called exponential-polynomial families (EPFs) [29]. PEFs ${\mathcal{E}}_D$ are regular minimal exponential families with polynomial sufficient statistics $t_i\left(x\right)=x^i$ for $i\in \{1,\dots ,D\}$. For example, the exponential distributions $\{p_{\lambda }\left(x\right)=\lambda exp(-\lambda x)\}$ form a PEF with $D=1$, $t\left(x\right)=x$ and $\mathcal{X}={\mathbb{R}}_{+}$, and the normal distributions form an EPF with $D=2$, $t\left(x\right)={\left[x{x}^2\right]}^{\top }$ and $\mathcal{X}=\mathbb{R}$, etc. Although the log-normalizer $F\left(\theta \right)$ can be obtained in closed-form for lower order PEFs (e.g., $D=1$ or $D=2$) or very special subfamilies (e.g., when $D=1$ and $t_1\left(x\right)=x^k$, exponential-monomial families [30]), a no-closed form formula is available for $F\left(\theta \right)$ of EPFs in general as soon $D\ge 4$ [31,32], and the cumulant function $F\left(\theta \right)$ is said to be computationally intractable. Notice that when $\mathcal{X}=\mathbb{R}$, the leading coefficient ${\theta }_D$ is negative for even integer order D. EPFs are attractive because these families can universally model any smooth multimodal distribution [28] and require fewer parameters in comparison to GMMs: Indeed, a univariate k-GMM $m\left(x\right)$ (at most k modes and $k-1$ antimodes) requires $3k-1$ parameters to specify $m\left(x\right)$ (or $k+1$ for a KDE with constant kernel width $\sigma$ or $2k-1$ for a KDE with varying kernel widths, but then $k=n$ observations). A density of an EPF of order D is called an exponential-polynomial density (EPD) and requires D parameters to specify $\theta$, with, at most, $\frac{D}{2}$ modes (and $\frac{D}{2}-1$ antimodes). The case of the quartic (polynomial) exponential densities ${\mathcal{E}}_4$ ($D=4$) has been extensively investigated in [31,33,34,35,36,37]. Armstrong and Brigo [38] discussed order-6 PEDs, and Efron and Hastie reported and order-7 PEF in their textbook (see Figure 5.7 of [39]). Figure 1 displays two examples of converting a GMM into a pair of dually parameterized exponential-polynomial densities.

Then by converting both mixture m and mixture $m^{\prime }$ into pairs of dually natural/moment parameterized unnormalized PEDs, i.e., $m\to (q_{{\overline{\theta }}_{\mathrm{SME}}},q_{{\overline{\eta }}_{\mathrm{MLE}}})$ and $m^{\prime }\to (q_{{\overline{\theta }}_{\mathrm{SME}}^{\prime }},q_{{\overline{\eta }}_{\mathrm{MLE}}}^{\prime })$, we approximate the JD between mixtures m and $m^{\prime }$ by using the four parameters of the PEDs

$$D_J[m,m^{\prime }]\approx {({\overline{\theta }}_{\mathrm{SME}}^{\prime }-{\overline{\theta }}_{\mathrm{SME}})}^{\top }({\overline{\eta }}_{\mathrm{MLE}}^{\prime }-{\overline{\eta }}_{\mathrm{MLE}}).$$

(8)

Let ${\mathrm{\Delta }}_J$ denote the approximation formula obtained from the two pairs of PEDs:

$${\mathrm{\Delta }}_J[p_{{\theta }_{\mathrm{SME}}},p^{{\eta }_{\mathrm{MLE}}};p_{{\theta }_{\mathrm{SME}}^{\prime }},p^{{\eta }_{\mathrm{MLE}}^{\prime }}]:={({\theta }_{\mathrm{SME}}^{\prime }-{\theta }_{\mathrm{SME}})}^{\top }({\eta }_{\mathrm{MLE}}^{\prime }-{\eta }_{\mathrm{MLE}}).$$

(9)

Let ${\mathrm{\Delta }}_J({\theta }_{\mathrm{SME}},{\eta }_{\mathrm{MLE}};{\theta }_{\mathrm{SME}}^{\prime },{\eta }_{\mathrm{MLE}}^{\prime }):={\mathrm{\Delta }}_J[p_{{\theta }_{\mathrm{SME}}},p^{{\eta }_{\mathrm{MLE}}};p_{{\theta }_{\mathrm{SME}}^{\prime }},p^{{\eta }_{\mathrm{MLE}}^{\prime }}]$. Then we have

$$D_J[m,m^{\prime }]\approx {\tilde{D}}_J[m,m^{\prime }]:={\mathrm{\Delta }}_J({\theta }_{\mathrm{SME}},{\eta }_{\mathrm{MLE}};{\theta }_{\mathrm{SME}}^{\prime },{\eta }_{\mathrm{MLE}}^{\prime }).$$

Note that ${\mathrm{\Delta }}_J$ is not a proper divergence as it may be negative since, in general, ${\overline{\eta }}_{\mathrm{MLE}}\ne \nabla F\left({\overline{\theta }}_{\mathrm{SME}}\right)$. That is, ${\mathrm{\Delta }}_J$ may not satisfy the law of the indiscernibles. Approximation ${\mathrm{\Delta }}_J$ is exact when $k_1=k_2=1$, with both m and $m^{\prime }$ belonging to an exponential family.

We experimentally show in Section 4 that the ${\tilde{D}}_J$ heuristic yields fast approximations of the JD compared to the MC baseline estimations by several order of magnitudes while approximating the JD reasonably well when the mixtures have a small number of modes.

For example, Figure 2 displays the unnormalized PEDs obtained for two Gaussian mixture models ($k_1=10$ components and $k_2=11$ components) into PEDs of a PEF of order $D=8$. The MC estimation of the JD with $s={10}^6$ samples yields $0.2633\cdots$, while the PED approximation of Equation (8) on corresponding PEFs yields $0.2618\dots$ (the relative error is $0.00585\cdots$ or about $0.585\dots \%$). It took about $2642.581$ milliseconds (with $s={10}^6$ on a Dell Inspiron 7472 laptop) to MC estimate the JD, while it took about $0.827$ milliseconds with the PEF approximation. Thus, we obtained a speed-up factor of about 3190 (three orders of magnitude) for this particular example. Notice that when viewing Figure 2, we tend to visually evaluate the dissimilarity using the total variation distance (a metric distance):

$$D_{\mathrm{TV}}[m,m^{\prime }]:=\frac{1}{2}\int |m\left(x\right)-m^{\prime }\left(x\right)|\mathrm{d}x,$$

rather than by a dissimilarity relating to the KLD. Using Pinsker’s inequality [40,41], we have $D_J[m,m^{\prime }]\ge D_{\mathrm{TV}}{[m,m^{\prime }]}^2$ and $D_{\mathrm{TV}}[m,m^{\prime }]\in [0,1]$. Thus, large TV distance (e.g., $D_{\mathrm{TV}}[m,m^{\prime }]=0.1$) between mixtures may have a small JD since Pinsker’s inequality yields $D_J[m,m^{\prime }]\ge 0.01$.

Let us point out that our approximation heuristic is deterministic, while the MC estimations are stochastic: That is, each MC run (Equation (4)) returns a different result, and a single MC run may yield a very bad approximation of the true Jeffreys divergence.

We compare our fast heuristic ${\tilde{D}}_J[m,m^{\prime }]={({\theta }_{\mathrm{SME}}^{\prime }-{\theta }_{\mathrm{SME}})}^{\top }({\eta }_{\mathrm{MLE}}^{\prime }-{\eta }_{\mathrm{MLE}})$ with two more costly methods relying on numerical procedures to convert natural ↔ moment parameters:

• Simplify GMMs $m_i$ into $p^{{\eta }_i^{\mathrm{MLE}}}$, and approximately convert the ${\overline{\eta }}_i^{\mathrm{MLE}}$’s into ${\tilde{\theta }}_i^{\mathrm{MLE}}$’s. Then approximate the Jeffreys divergence as

  $$D_J[m_1,m_2]\simeq {\tilde{\mathrm{\Delta }}}_J^{\mathrm{MLE}}[m_1,m_2]:={({\tilde{\theta }}_2^{\mathrm{MLE}}-{\tilde{\theta }}_1^{\mathrm{MLE}})}^{\top }({\overline{\eta }}_2^{\mathrm{MLE}}-{\overline{\eta }}_1^{\mathrm{MLE}}).$$

  (10)
• Simplify GMMs $m_i$ into $p_{{\overline{\theta }}_i^{\mathrm{SME}}}$, and approximately convert the ${\overline{\theta }}_i^{\mathrm{SME}}$’s into ${\tilde{\eta }}_i^{\mathrm{SME}}$’s. Then approximate the Jeffreys divergence as

  $$D_J[m_1,m_2]\simeq {\tilde{\mathrm{\Delta }}}_J^{\mathrm{SME}}(m_1,m_2)={({\overline{\theta }}_2^{\mathrm{SME}}-{\overline{\theta }}_1^{\mathrm{SME}})}^{\top }({\tilde{\eta }}_2^{\mathrm{SME}}-{\tilde{\eta }}_1^{\mathrm{SME}}).$$

  (11)

1.4. Contributions and Paper Outline

Our contributions are summarized as follows:

• We explain how to convert any continuous density $r\left(x\right)$ (including GMMs) into a polynomial exponential density in Section 2 using integral-based extensions of the Maximum Likelihood Estimator [22] (MLE estimates in the moment parameter space H, Theorem 1 and Corollary 1) and the Score Matching Estimator [27] (SME estimates in the natural parameter space $\Theta$, Theorem 3). We show a connection between SME and the Moment Linear System Estimator [28] (MLSE).
• We report a closed-form formula to evaluate the goodness-of-fit of a polynomial family density to a GMM in Section 3 using an extension of the Hyvärinen divergence [42] (Theorem 4) and discuss the problem of model selection for choosing the order D of the polynomial exponential family.
• We show how to approximate the Jeffreys divergence between GMMs using a pair of natural/moment parameter PED conversion and present experimental results that display a gain of several orders of magnitude of performance when compared to the vanilla Monte Carlo estimator in Section 4. We observe that the quality of the approximations depend on the number of modes of the GMMs [43]. However, calculating or counting the modes of a GMM is a difficult problem in its own [43].

The paper is organized as follows: In Section 2, we show how to convert arbitrary probability density functions into polynomial exponential densities using the integral-based Maximum Likelihood Estimator (MLE) and Score Matching Estimator (SME). We describe a Maximum Entropy method to iteratively convert moment parameters into natural parameters in Section 2.3.1. It is followed by Section 3, which shows how to calculate in closed-form the order-2 Hyvärinen divergence between a GMM and a polynomial exponential density. We use this criterion to perform model selection. Section 4 presents our computational experiments that demonstrate a gain of several orders of magnitudes for GMMs with a small number of modes. Finally, we conclude in Section 5.

2. Converting Finite Mixtures to Exponential Family Densities

We report two generic methods to convert a mixture $m\left(x\right)$ into a density $p_{\theta }\left(x\right)$ of an exponential family: The first method extending the MLE in Section 2.1 proceeds using the mean parameterization $\eta$, while the second method extending the SME in Section 2.2 uses the natural parameterization of the exponential family. We then describe how to convert the moments parameters into natural parameters (and vice versa) for polynomial exponential families in Section 2.3. We show how to instantiate these generic conversion methods for GMMs: It requires calculating non-central moments of GMMs in closed-form. The efficient computations of raw moments of GMMs is detailed in Section 2.4.

2.1. Conversion Using the Moment Parameterization (MLE)

Let us recall that in order to estimate the moment or mean parameter ${\widehat{\eta }}_{\mathrm{MLE}}$ of a density belonging an exponential family

$${\mathcal{E}}_t:=\left\{p_{\theta }\left(x\right)=exp\left(t{\left(x\right)}^{\top }\theta -F\left(\theta \right)\right)\right\}$$

with a sufficient statistic vector $t\left(x\right)={[t_1\left(x\right)\dots t_D\left(x\right)]}^{\top }$ from an i.i.d. sample set $x_1,\dots ,x_n$, the Maximum Likelihood Estimator (MLE) [22,44] yields

$$\begin{array}{ccc}& & \underset{\theta }{max}\prod _{i=1}^n{p}_{\theta }\left(x_i\right),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& \equiv & \underset{\theta }{max}\sum _{i=1}^{n}log{p}_{\theta }\left(x_i\right),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& =& \underset{\theta }{max}E\left(\theta \right):=\left(\sum _{i=1}^{n}t{\left(x_i\right)}^{\top }\theta \right)-nF\left(\theta \right),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& \Rightarrow & {\widehat{\eta }}_{\mathrm{MLE}}=\frac{1}{n}\sum _{i=1}^{n}t\left(x_i\right).\hfill \end{array}$$

(15)

In statistics, Equation (12) is called the estimating equation. The MLE exists under mild conditions [22] and is unique since the Hessian ${\nabla }^{2}E\left(\theta \right)={\nabla }^{2}F\left(\theta \right)$ of the estimating equation is positive-definite (log-normalizers $F\left(\theta \right)$ are always strictly convex and real analytic [22]). The MLE is consistent and asymptotically normally distributed [22]. Furthermore, since the MLE satisfies the equivariance property [22], we have ${\widehat{\theta }}_{\mathrm{MLE}}=\nabla F^{*}\left({\widehat{\eta }}_{\mathrm{MLE}}\right)$, where $\nabla F^{*}$ denotes the gradient of the conjugate function $F^{*}\left(\eta \right)$ of the cumulant function $F\left(\theta \right)$ of the exponential family. In general, $\nabla F^{*}$ is intractable for PEDs with $D\ge 4$.

By considering the empirical distribution

$$p_e\left(x\right):=\frac{1}{n}\sum _{i=1}^s{\delta }_{x_i}\left(x\right),$$

where ${\delta }_{x_i}(·)$ denotes the Dirac distribution at location $x_i$, we can formulate the MLE problem as a minimum KLD problem between the empirical distribution and a density of the exponential family:

$$\begin{array}{ccc}\hfill \underset{\theta }{min}{D}_{\mathrm{KL}}[p_e:p_{\theta }]& =& min-H\left[p_e\right]-E_{p_e}[log{p}_{\theta }\left(x\right)],\hfill \\ & \equiv & \underset{\theta }{max}\frac{1}{n}\sum _{i=1}^{n}log{p}_{\theta }\left(x_i\right),\hfill \end{array}$$

since the entropy term $H\left[p_e\right]$ is independent of $\theta$.

Thus, to convert an arbitrary smooth density $r\left(x\right)$ into a density $p_{\theta }$ of an exponential family ${\mathcal{E}}_t$, we have to solve the following minimization problem:

$$\underset{\theta \in \Theta }{min}{D}_{\mathrm{KL}}[r:p_{\theta }].$$

Rewriting the minimization problem as:

$$\begin{array}{ccc}& & \underset{\theta }{min}{D}_{\mathrm{KL}}[r:p_{\theta }]=-\int r\left(x\right)log{p}_{\theta }\left(x\right)\mathrm{d}x+\int r\left(x\right)logr\left(x\right)\mathrm{d}x,\hfill \\ & \equiv & \underset{\theta }{min}-\int r\left(x\right)log{p}_{\theta }\left(x\right)\mathrm{d}x,\hfill \\ & =& \underset{\theta }{min}\int r\left(x\right)(F\left(\theta \right)-{\theta }^{\top }t\left(x\right))\mathrm{d}x,\hfill \\ & =& \underset{\theta }{min}\overline{E}\left(\theta \right)=F\left(\theta \right)-{\theta }^{\top }{E}_r\left[t\left(x\right)\right],\hfill \end{array}$$

we obtain

$${\overline{\eta }}_{\mathrm{MLE}}\left(r\right):=E_r\left[t\left(x\right)\right]={\int }_{\mathcal{X}}r\left(x\right)t\left(x\right)\mathrm{d}x.$$

(16)

The minimum is unique since ${\nabla }^2\overline{E}\left(\theta \right)={\nabla }^{2}F\left(\theta \right)\succ 0$ (positive-definite matrix). This conversion procedure $r\left(x\right)\to p^{{\overline{\eta }}_{\mathrm{MLE}}\left(r\right)}\left(x\right)$ can be interpreted as an integral extension of the MLE, hence the $\overline{\dot{}}$ notation in ${\overline{\eta }}_{\mathrm{MLE}}$. Notice that the ordinary MLE is ${\widehat{\eta }}_{\mathrm{MLE}}={\overline{\eta }}_{\mathrm{MLE}}\left(p_e\right)$ obtained for the empirical distribution: $r=p_e$: ${\overline{\eta }}_{\mathrm{MLE}}\left(p_e\right)=\frac{1}{n}{\sum }_{i=1}^{n}t\left(x_i\right)$.

Theorem 1.

The best density $p^{\overline{\eta }}\left(x\right)$ of an exponential family ${\mathcal{E}}_t=\{p_{\theta }:\theta \in \Theta \}$ minimizing the Kullback–Leibler divergence $D_{\mathrm{KL}}[r:p_{\theta }]$ between a density r and a density $p_{\theta }$ of an exponential family ${\mathcal{E}}_t$ is $\overline{\eta }=E_r\left[t\left(x\right)\right]={\int }_{\mathcal{X}}r\left(x\right)t\left(x\right)\mathrm{d}x$.

Notice that when $r=p_{\theta }$, we obtain $\overline{\eta }=E_{p_{\theta }}\left[t\left(x\right)\right]=\eta$ so that the method ${\overline{\eta }}_{\mathrm{MLE}}\left(r\right)$ is consistent (by analogy to the finite i.i.d. MLE case): ${\overline{\eta }}_{\mathrm{MLE}}\left(p_{\theta }\right)=\eta =\nabla F\left(\theta \right)$.

The KLD right-sided minimization problem can be interpreted as an information projection of r onto ${\mathcal{E}}_t$. As a corollary of Theorem 1, we obtain:

Corollary 1

(Best right-sided KLD simplification of a mixture). The best right-sided KLD simplification of a homogeneous mixture of exponential families [2] $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_{{\theta }_i}\left(x\right)$ with $p_{{\theta }_i}\in {\mathcal{E}}_t$, i.e., ${min}_{\theta \in \Theta }{D}_{\mathrm{KL}}[m:p_{\theta }]$, into a single component $p^{\eta }\left(x\right)$ is given by $\eta ={\widehat{\eta }}_{\mathrm{MLE}}\left(m\right)=E_m\left[t\left(x\right)\right]={\sum }_{i=1}^k{\eta }_i=\overline{\eta }$.

Equation (16) allows us to greatly simplify the proofs reported in [45] for mixture simplifications that involved the explicit use of the Pythagoras’ theorem in the dually flat spaces of exponential families [42]. Figure 3 displays the geometric interpretation of the best KLD simplification of a GMM with ambient space the probability space $(\mathbb{R},\mathcal{B}\left(\mathbb{R}\right),{\mu }_L)$, where ${\mu }_L$ denotes the Lebesgue measure and $\mathcal{B}\left(\mathbb{R}\right)$ the Borel $\sigma$-algebra of $\mathbb{R}$.

Let us notice that Theorem 1 yields an algebraic system for polynomial exponential densities, i.e., $E_m\left[x^i\right]={\overline{\eta }}_i$ for $i\in \{1,\dots ,D\}$, to compute ${\overline{\eta }}_{\mathrm{MLE}}\left(m\right)$ for a given GMM $m\left(x\right)$ (since raw moments $E_m\left[x^i\right]$ are algebraic). In contrast with this result, the MLE of i.i.d. observations is in general not an algebraic function [46] but a transcendental function.

2.2. Converting to a PEF Using the Natural Parameterization (SME)

Integral-Based Score Matching Estimator (SME)

To convert the density $r\left(x\right)$ into an exponential density with sufficient statistics $t\left(x\right)$, we can also use the Score Matching Estimator [27,47] (SME). The Score Matching Estimator minimizes the Hyvärinen divergence $D_H$ (Equation (4) of [47]):

$$D_H[p:p_{\theta }]:=\frac{1}{2}\int {\parallel {\nabla }_{x}logp\left(x\right)-{\nabla }_{x}log{p}_{\theta }\left(x\right)\parallel }^{2}p\left(x\right)\mathrm{d}x.$$

The Hyvärinen divergence is also known as half of the relative Fisher information in the optimal transport community (Equation (8) of [48] or Equation (2.2) in [49]), where it is defined for two measures $\mu$ and $\nu$ as follows:

$$I[\mu :\nu ]:={\int }_{\mathcal{X}}{∥\nabla log\frac{\mathrm{d}\mu }{\mathrm{d}\nu }∥}^2\mathrm{d}\mu =4{\int }_{\mathcal{X}}{∥\nabla \sqrt{\frac{\mathrm{d}\mu }{\mathrm{d}\nu }}∥}^2\mathrm{d}\nu .$$

Moreover, the relative Fisher information can be defined on complete Riemannian manifolds [48].

That is, we convert a density $r\left(x\right)$ into an exponential family density $p_{\theta }\left(x\right)$ using the following minimizing problem:

$${\theta }_{\mathrm{SME}}\left(r\right)=\underset{\theta \in \Theta }{min}{D}_H[r:p_{\theta }].$$

Beware that in statistics, the score $s_{\theta }\left(x\right)$ is defined by ${\nabla }_{\theta }log{p}_{\theta }\left(x\right)$, but in Score Matching, we refer to the “data score” defined by ${\nabla }_{x}log{p}_{\theta }\left(x\right)$. Hyvärinen [47] gave an explanation of the naming “score” using a spurious location parameter.

• Generic solution: It can be shown that for exponential families [47], we obtain the following solution:

  $${\theta }_{\mathrm{SME}}\left(r\right)=-{\left(E_r\left[A\left(x\right)\right]\right)}^{-1}\times \left(E_r\left[b\left(x\right)\right]\right),$$

  (17)

  where

  $$A\left(x\right):={\left[t_i^{\prime }\left(x\right)t_j^{\prime }\left(x\right)\right]}_{ij}$$

  is a $D\times D$ symmetric matrix, and

  $$b\left(x\right)={[t_1^{″}\left(x\right)\dots t_D^{″}\left(x\right)]}^{\top }$$

  is a D-dimensional column vector.

Theorem 2. 

The best conversion of a density $r\left(x\right)$ into a density $p_{\theta }\left(x\right)$ of an exponential family minimizing the right-sided Hyvärinen divergence is

$${\theta }_{\mathrm{SME}}\left(r\right)=-{\left(E_r\left[{\left[t_i^{\prime }\left(x\right)t_j^{\prime }\left(x\right)\right]}_{ij}\right]\right)}^{-1}\times \left(E_r{\left[[t_1^{″}\left(x\right)\dots t_D^{″}\left(x\right)]\right]}^{\top }\right).$$

• Solution instantiated for polynomial exponential families:
  For polynomial exponential families of order D, we have $t_i^{\prime }\left(x\right)=i{x}^{i-1}$ and $t_i^{″}\left(x\right)=i(i-1)x^{i-2}$, and therefore, we have

  $$A_D=E_r\left[A\left(x\right)\right]={\left[ij{\mu }_{i+j-2}\left(r\right)\right]}_{ij},$$

  and

  $$b_D=E_s\left[b\left(x\right)\right]={\left[j(j-1){\mu }_{j-2}\left(r\right)\right]}_j,$$

  where ${\mu }_l\left(r\right):=E_r\left[X^l\right]$ denotes the l-th raw moment of distribution $X\sim r\left(x\right)$ (with the convention that $m_{-1}\left(r\right)=0$). For a probability density function $r\left(x\right)$, we have ${\mu }_1\left(r\right)=1$.
  Thus, the integral-based SME of a density r is:

  $${\theta }_{\mathrm{SME}}\left(r\right)=-{\left({\left[ij{\mu }_{i+j-2}\left(r\right)\right]}_{ij}\right)}^{-1}\times {\left[j(j-1){\mu }_{j-2}\left(r\right)\right]}_j.$$

  (18)
  For example, matrix $A_4$ is

  $$\left[\begin{array}{cccc}{\mu }_0& 2{\mu }_1& 3{\mu }_2& 4{\mu }_3\\ 2{\mu }_1& 4{\mu }_2& 6{\mu }_3& 8{\mu }_4\\ 3{\mu }_2& 6{\mu }_3& 9{\mu }_4& 12{\mu }_5\\ 4{\mu }_3& 8{\mu }_4& 12{\mu }_5& 16{\mu }_6\end{array}\right].$$
• Faster PEF solutions using Hankel matrices:
  The method of Cobb et al. [28] (1983) anticipated the Score Matching method of Hyvärinen (2005). It can be derived from Stein’s lemma for exponential families [50]. The integral-based Score Matching method is consistent, i.e., if $r=p_{\theta }$, then ${\overline{\theta }}_{\mathrm{SME}}=\theta$: The probabilistic proof for $r\left(x\right)=p_e\left(x\right)$ is reported as Theorem 2 of [28]. The integral-based proof is based on the property that arbitrary order partial mixed derivatives can be obtained from higher-order partial derivatives with respect to ${\theta }_1$ [29]:

  $${\partial }_1^{i_1}\dots {\partial }_D^{i_D}F\left(\theta \right)={\partial }_1^{{\sum }_{j=1}^{D}j{i}_j}F\left(\theta \right),$$

  where ${\partial }_i:=\frac{\partial }{\partial {\theta }_i}$.
  The complexity of the direct SME method is $O\left(D^3\right)$ as it requires the inverse of the $D\times D$-dimensional matrix $A_D$.
  We show how to lower this complexity by reporting an equivalent method (originally presented in [28]) that relies on recurrence relationships between the moments of $p_{\theta }\left(x\right)$ for PEDs. Recall that ${\mu }_l\left(r\right)$ denotes the l-th raw moment $E_r\left[x^l\right]$.
  Let $A^{\prime }={\left[a_{i+j-2}^{\prime }\right]}_{ij}$ denote the $D\times D$ symmetric matrix with $a_{i+j-2}^{\prime }\left(r\right)={\mu }_{i+j-2}\left(r\right)$ (with $a_0^{\prime }\left(r\right)={\mu }_0\left(r\right)=1$), and $b^{\prime }={\left[b_i\right]}_i$ the D-dimensional vector with $b_i^{\prime }\left(r\right)=(i+1){\mu }_i\left(r\right)$. We solve the system $A^{\prime }\beta =b^{\prime }$ to obtain $\beta ={A^{\prime }}^{-1}{b}^{\prime }$. We then obtain the natural parameter ${\overline{\theta }}_{\mathrm{SME}}$ from the vector $\beta$ as

  $${\overline{\theta }}_{\mathrm{SME}}=\left[\begin{array}{c}-\frac{{\beta }_1}{2}\\ ⋮\\ -\frac{{\beta }_i}{i+1}\\ ⋮\\ -\frac{{\beta }_D}{D+1}\end{array}\right].$$

  (19)
  Now, if we inspect matrix $A_D^{\prime }=\left[{\mu }_{i+j-2}\left(r\right)\right]$, we find that matrix $A_D^{\prime }$ is a Hankel matrix: A Hankel matrix has constant anti-diagonals and can be inverted in quadratic-time [51,52] instead of cubic time for a general $D\times D$ matrix. (The inverse of a Hankel matrix is a Bezoutian matrix [53].) Moreover, a Hankel matrix can be stored using linear memory (store $2D-1$ coefficients) instead of quadratic memory of regular matrices.
  For example, matrix $A_4^{\prime }$ is:

  $$A_4^{\prime }=\left[\begin{array}{cccc}{\mu }_0& {\mu }_1& {\mu }_2& {\mu }_3\\ {\mu }_1& {\mu }_2& {\mu }_3& {\mu }_4\\ {\mu }_2& {\mu }_3& {\mu }_4& {\mu }_5\\ {\mu }_3& {\mu }_4& {\mu }_5& {\mu }_6\end{array}\right],$$

  and requires only $6=2\times 4-2$ coefficients to be stored instead of $4\times 4=16$. The order-d moment matrix is

  $$A_d^{\prime }:={\left[{\mu }_{i+j-2}\right]}_{ij}=\left[\begin{array}{cccc}{\mu }_0& {\mu }_1& \dots & {\mu }_d\\ {\mu }_1& {\mu }_2& \dots & ⋮\\ ⋮& & \ddots & ⋮\\ {\mu }_d& \dots & \dots & {\mu }_{2d}\end{array}\right],$$

  is a Hankel matrix stored using $2d+1$ coefficients:

  $$A_d^{\prime }=:\mathrm{Hankel}({\mu }_0,{\mu }_1,\dots ,{\mu }_{2d}).$$
  In statistics, those matrices $A_d^{\prime }$ are called moment matrices and well-studied [54,55,56]. The variance $\mathrm{Var}\left[X\right]$ of a random variable X can be expressed as the determinant of the order-2 moment matrix:

  $$\mathrm{Var}\left[X\right]=E\left[{(X-\mu )}^2\right]=E\left[X^2\right]-E{\left[X\right]}^2={\mu }_2-{\mu }_1^2=\det \left(\left[\begin{array}{cc}1& {\mu }_1\\ {\mu }_1& {\mu }_2\end{array}\right]\right)\ge 0.$$
  This observation yields a generalization of the notion of variance to $d+1$ random variables: $X_1,\dots ,X_{d+1}{\sim }_{iid}{F}_X\Rightarrow E\left[{\prod }_{j>i}{(X_i-X_j)}^2\right]=(d+1)!\det \left(M_d\right)\ge 0$. The variance can be expressed as $E\left[\frac{1}{2}{(X_1-X_2)}^2\right]$ for $X_1,X_2{\sim }_{iid}{F}_X$. See [57] (Chapter 5) for a detailed description related to U-statistics.
  For GMMs r, the raw moments ${\mu }_l\left(r\right)$ to build matrix $A_D$ can be calculated in closed-form, as explained in Section 2.4.

Theorem 3 (Score matching GMM conversion)

The Score Matching conversion of a GMM $m\left(x\right)$ into a polynomial exponential density $p_{{\theta }_{\mathrm{SME}}\left(m\right)}\left(x\right)$ of order D is obtained as

$${\theta }_{\mathrm{SME}}\left(m\right)=-{\left({\left[ij{m}_{i+j-2}\right]}_{ij}\right)}^{-1}\times {\left[j(j-1)m_{j-2}\right]}_j,$$

where $m_i=E_m\left[x^i\right]$ denote the ith non-central moment of the GMM $m\left(x\right)$.

2.3. Converting Numerically Moment Parameters from/to Natural Parameters

Recall that our fast heuristic approximates the Jeffreys divergence by

$${\tilde{D}}_J[m,m^{\prime }]:={({\overline{\theta }}_{\mathrm{SME}}\left(m^{\prime }\right)-{\overline{\theta }}_{\mathrm{SME}}\left(m\right))}^{\top }({\overline{\eta }}_{\mathrm{MLE}}\left(m^{\prime }\right)-{\overline{\eta }}_{\mathrm{MLE}}\left(m\right)).$$

Because F and $\nabla F^{*}$ are not available in closed form (except for the case $D=2$ of the normal family), we cannot obtain $\theta$ from a given $\eta$ (using $\theta =\nabla F^{*}\left(\eta \right)$) nor $\eta$ from a given $\theta$ (using $\eta =\nabla F\left(\theta \right)$).

However, provided that we can approximate numerically $\tilde{\eta }\simeq \nabla F\left(\theta \right)$ and $\tilde{\theta }\simeq \nabla F^{*}\left(\eta \right)$, we also consider these two approximations for the Jeffreys divergence:

$${\tilde{\mathrm{\Delta }}}_J^{\mathrm{MLE}}[m_1,m_2]:={({\tilde{\theta }}_2^{\mathrm{MLE}}-{\tilde{\theta }}_1^{\mathrm{MLE}})}^{\top }({\overline{\eta }}_2^{\mathrm{MLE}}-{\overline{\eta }}_1^{\mathrm{MLE}}),$$

and

$${\tilde{\mathrm{\Delta }}}_J^{\mathrm{SME}}[m_1,m_2]={({\overline{\theta }}_2^{\mathrm{SME}}-{\overline{\theta }}_1^{\mathrm{SME}})}^{\top }({\tilde{\eta }}_2^{\mathrm{SME}}-{\tilde{\eta }}_1^{\mathrm{SME}}).$$

We show how to numerically estimate ${\tilde{\theta }}^{\mathrm{MLE}}\simeq \nabla F\left({\overline{\eta }}^{\mathrm{MLE}}\right)$ from ${\overline{\eta }}^{\mathrm{MLE}}$ in Section 2.3.1. Next, in Section 2.3.2, we show how to stochastically estimate ${\tilde{\eta }}^{\mathrm{SME}}\simeq \nabla F^{*}\left({\overline{\theta }}^{\mathrm{SME}}\right)$.

2.3.1. Converting Moment Parameters to Natural Parameters Using Maximum Entropy

Let us report the iterative approximation technique of [58] (which extended the method described in [35]) based on solving a maximum entropy problem (MaxEnt problem). This method will be useful when comparing our fast heuristic ${\tilde{D}}_J[m,m^{\prime }]$ with the approximations ${\tilde{\mathrm{\Delta }}}_J^{\mathrm{MLE}}[m,m^{\prime }]$ and ${\tilde{\mathrm{\Delta }}}_J^{\mathrm{SME}}[m,m^{\prime }]$.

The density $p_{\theta }$ of any exponential family can be characterized as a maximum entropy distribution given the D moment constraints $E_{p_{\theta }}\left[t_i\left(x\right)\right]={\eta }_i$: Namely, ${max}_{p}h\left(p\right)$ subject to the $D+1$ moment constraints $\int t_i\left(x\right)p\left(x\right)\mathrm{d}x={\eta }_i$ for $i\in \{0,\dots ,D\}$, where we added by convention ${\eta }_0=1$ and $t_0\left(x\right)=1$ (so that $\int p\left(x\right)\mathrm{d}x=1$). The solution of this MaxEnt problem [58] is $p\left(x\right)=p_{\lambda }$, where $\lambda$ are the $D+1$ Lagrangian parameters. Here, we adopt the following canonical parameterization of the densities of an exponential family:

$$p_{\lambda }\left(x\right):=exp\left(-\sum _{i=0}^D{\lambda }_i{t}_i\left(x\right)\right).$$

That is, $F\left(\lambda \right)={\lambda }_0$ and ${\lambda }_i=-{\theta }_i$ for $i\in \{1,\dots ,D\}$. Parameter $\lambda$ is a kind of augmented natural parameter that includes the log-normalizer in its first coefficient.

Let $K_i\left(\lambda \right):=E_{p_{\theta }}\left[t_i\left(x\right)\right]={\eta }_i$ denote the set of $D+1$ non-linear equations for $i\in \{0,\dots ,D\}$. The Iterative Linear System Method [58] (ILSM) converts $p^{\eta }$ to $p_{\theta }$ iteratively. We initialize ${\lambda }^{\left(0\right)}$ to ${\overline{\theta }}_{\mathrm{SME}}$ (and calculate numerically ${\lambda }_0^{\left(0\right)}=F\left({\overline{\theta }}_{\mathrm{SME}}\right)$).

At iteration t with current estimate ${\lambda }^{\left(t\right)}$, we use the following first-order Taylor approximation:

$$K_i\left(\lambda \right)\approx K_i\left({\lambda }^{\left(t\right)}\right)+(\lambda -{\lambda }^{\left(t\right)})\nabla K_i\left({\lambda }^{\left(t\right)}\right).$$

Let $H\left(\lambda \right)$ denote the $(D+1)\times (D+1)$ matrix:

$$H\left(\lambda \right):={\left[\frac{\partial K_i\left(\lambda \right)}{\partial {\theta }_j}\right]}_{ij}.$$

We have

$$H_{ij}\left(\lambda \right)=H_{ji}\left(\lambda \right)=-E_{p_{\theta }}\left[t_i\left(x\right)t_j\left(x\right)\right].$$

We update as follows:

$${\lambda }^{(t+1)}={\lambda }^{\left(t\right)}+H^{-1}\left({\lambda }^{\left(t\right)}\right)\left[\begin{array}{c}{\eta }_0-K_0\left({\lambda }^{\left(t\right)}\right)\\ ⋮\\ {\eta }_D-K_D\left({\lambda }^{\left(t\right)}\right)\end{array}\right].$$

(20)

For a PEF of order D, we have

$$H_{ij}\left(\lambda \right)=-E_{p_{\theta }}\left[x^{i+j-2}\right]=-{\mu }_{i+j-2}\left(p_{\theta }\right).$$

This yields a moment matrix $H_{\lambda }$ (Hankel matrix), which can be inverted in quadratic time [52]. In our setting, the moment matrix is invertible because $\left|H\right|>0$, see [59].

Let ${\tilde{\lambda }}_T\left(\eta \right)$ denote ${\theta }^{\left(T\right)}$ after T iterations (retrieved from ${\lambda }^{\left(T\right)}$) and the corresponding natural parameter of the PED. We have the following approximation of the JD:

$$D_J[m,m^{\prime }]\approx {({\tilde{\theta }}_T\left({\eta }^{\prime }\right)-{\tilde{\theta }}_T\left(\eta \right))}^{\top }({\eta }^{\prime }-\eta ).$$

The method is costly because we need to numerically calculate ${\mu }_{i+j-2}\left(p_{\theta }\right)$ and the $K_i$’s (e.g., univariate Simpson integrator). Another potential method consists of estimating these expectations using acceptance-rejection sampling [60,61]. We may also consider the holonomic gradient descent [29]. Thus, the conversion $\eta \to \theta$ method is costly. Our heuristic ${\tilde{\mathrm{\Delta }}}_J$ bypasses this costly moment-to-natural parameter conversion by converting each mixture m to a pair $(p_{{\theta }_{\mathrm{SME}}},p_{{\eta }_{\mathrm{MLE}}})$ of PEDs parameterized in the natural and moment parameters (i.e., loosely speaking, we untangle these dual parameterizations).

2.3.2. Converting Natural Parameters to Moment Parameters

Given a PED $p_{\theta }\left(x\right)$, we have to find its corresponding moment parameter $\eta$ (i.e., $p_{\theta }=p^{\eta }$). Since $\eta =E_{p_{\theta }}\left[t\left(x\right)\right]$, we sample s i.i.d. variates $x_1,\dots ,x_s$ from $p_{\theta }$ using acceptance-rejection sampling [60,61] or any other Markov chain Monte Carlo technique [62] and estimate $\widehat{\eta }$ as:

$$\widehat{\eta }=\frac{1}{s}\sum _{i=1}^{s}t\left(x_i\right).$$

2.4. Raw Non-Central Moments of Normal Distributions and GMMs

In order to implement the MLE or SME Gaussian mixture conversion procedures, we need to calculate the raw moments of a Gaussian mixture model. The l-th moment raw moment $E\left[Z^l\right]$ of a standard normal distribution $Z\sim N(0,1)$ is 0 when l is odd (since the normal standard density is an even function) and $(l-1)!!=2^{-\frac{l}{2}}\frac{l!}{(l/2)!}$ when l is even, where $n!!=\sqrt{\frac{2^{n+1}}{\pi }}\Gamma (\frac{n}{2}+1)={\prod }_{k=0}^{\lceil \frac{n}{2}\rceil -1}(n-2k)$ is the double factorial (with $(-1)!!=1$ by convention). Using the binomial theorem, we deduce that a normal distribution $X=\mu +\sigma Z$ has finite moments:

$${\mu }_l\left(p_{\mu ,\sigma }\right)=E_{p_{\mu ,\sigma }}\left[X^l\right]=E\left[{(\mu +\sigma Z)}^l\right]=E\left[{(\mu +\sigma Z)}^l\right]=\sum _{i=0}^l\left(\genfrac{}{}{0pt}{}{l}{i}\right){\mu }^{l-i}{\sigma }^{i}E\left[Z^i\right].$$

That is, we have

$${\mu }_l\left(p_{\mu ,\sigma }\right)=\sum _{i=0}^{\lfloor \frac{l}{2}\rfloor }\left(\genfrac{}{}{0pt}{}{l}{i}\right)(2i-1)!!{\mu }^{l-2i}{\sigma }^{2i},$$

(21)

where $n!!$ denotes the double factorial:

$$n!!=\prod _{k=0}^{\left.⌈\frac{n}{2}\right]-1}(n-2k)=\left\{\begin{array}{cc}\prod _{k=1}^{\frac{n}{2}}\left(2k\right)\hfill & n \mathrm{is} \mathrm{even},\hfill \\ \prod _{k=1}^{\frac{n+1}{2}}(2k-1)\hfill & n \mathrm{is} \mathrm{odd}.\hfill \end{array}\right.$$

By the linearity of the expectation $E[·]$, we deduce the l-th raw moment of a GMM $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_{{\mu }_i,{\sigma }_i}\left(x\right)$:

$${\mu }_l\left(m\right)=\sum _{i=1}^k{w}_i{\mu }_l\left(p_{{\mu }_I,{\sigma }_i}\right).$$

Notice that by using [63], we can extend this formula to truncated normals and GMMs. Thus, computing the first $O\left(D\right)$ raw moments of a GMM with k components can be done in $O\left(k{D}^2\right)$ using the Pascal triangle method for computing the binomial coefficients. See also [64].

3. Goodness-of-Fit between GMMs and PEDs: Higher Order Hyvärinen Divergences

Once we have converted a GMM $m\left(x\right)$ into an unnormalized PED $q_{{\theta }_m}\left(x\right)={\tilde{p}}_{{\theta }_m}\left(x\right)$, we would like to evaluate the quality of the conversion, i.e., $D[m\left(x\right):q_{{\theta }_m}\left(x\right)]$, using a statistical divergence $D[·:·]$. This divergence shall allow us to perform model selection by choosing the order D of the PEF so that $D[m\left(x\right):p_{\theta }\left(x\right)]\le ϵ$ for $\theta \in {\mathbb{R}}^D$, where $ϵ>0$ is a prescribed threshold. Since PEDs have computationally intractable normalization constants, we consider a right-sided projective divergence [42] $D[p:q]$ that satisfies $D[p:\lambda q]=D[p:q]=D[p:\tilde{q}]$ for any $\lambda >0$. For example, we may consider the $\gamma$-divergence [65] that is a two-sided projective divergence: $D_{\gamma }[\lambda p:{\lambda }^{\prime }q]=D[p:q]=D[\tilde{p}:\tilde{q}]$ for any $\lambda ,{\lambda }^{\prime }>0$ and converge to the KLD when $\gamma \to 0$. However, the $\gamma$-divergence between a mixture model and an unnormalized PEF does not yield a closed-form formula. Moreover, the $\gamma$-divergence between two unnormalized PEDs is expressed using the log-normalizer function $F(·)$ that is computationally intractable [66].

In order to a get a closed-form formula for a divergence between a mixture model and an unnormalized PED, we consider the order-$\alpha$ (for $\alpha >0$) Hyvärinen divergence [42] as follows:

$$D_{H,\alpha }[p:q]:=\int p{\left(x\right)}^{\alpha }{\left({\nabla }_{x}logp\left(x\right)-{\nabla }_{x}logq\left(x\right)\right)}^2\mathrm{d}x,\alpha >0.$$

(22)

The Hyvärinen divergence [42] (order-1 Hyvärinen divergence) has also been called the Fisher divergence [27,67,68,69] or relative Fisher information [48]. Notice that when $\alpha =1$, $D_{H,1}[p:q]=D_H[p:q]$, the ordinary Hyvärinen divergence [27].

The Hyvärinen divergences $D_{H,\alpha }$ is a right-sided projective divergence, meaning that the divergence satisfies $D_{H,\alpha }[p:q]=D_{H,\alpha }[p:\lambda q]$ for any $\lambda >0$. That is, we have $D_{H,\alpha }[p:q]=D_{H,\alpha }[p:\tilde{q}]$. Thus, we have $D_{H,\alpha }[m:p_{\theta }]=D_{H,\alpha }[m:q_{\theta }]$ for an unnormalized PED $q_{\theta }={\tilde{p}}_{\theta }$. For statistical estimation, it is enough to have a sided projective divergence since we need to evaluate the goodness of fit between the (normalized) empirical distribution $p_e$ and the (unnormalized) parameteric density.

For univariate distributions, ${\nabla }_{x}logp\left(x\right)=\frac{p^{\prime }\left(x\right)}{p\left(x\right)}$, and $\frac{p^{\prime }\left(x\right)}{p\left(x\right)}=\frac{{\tilde{p}}^{\prime }\left(x\right)}{\tilde{p}\left(x\right)}$, where $\tilde{p}\left(x\right)$ is the unnormalized model.

Let $P_{\theta }\left(x\right):={\sum }_{i=1}^D{\theta }_i{x}^i$ be a homogeneous polynomial defining the shape of the EPF:

$$p_{\theta }\left(x\right)=exp\left(P_{\theta }\left(x\right)-F\left(\theta \right)\right).$$

For PEDs with the homogeneous polynomial $P_{\theta }\left(x\right)$, we have $\frac{p^{\prime }\left(x\right)}{p\left(x\right)}={(log{P}_{\theta }\left(x\right))}^{\prime }={\sum }_{i=1}^{D}i{\theta }_i{x}^{i-1}$.

Theorem 4.

The Hyvärinen divergence $D_{H,2}[m:q_{\theta }]$ of order 2 between a Gaussian mixture $m\left(x\right)$ and a polynomial exponential family density $q_{\theta }\left(x\right)$ is available in closed form.

Proof. 

We have $D_{H,2}[m:q]=\int m{\left(x\right)}^2{\left(\frac{m^{\prime }\left(x\right)}{m\left(x\right)}-{\sum }_{i=1}^{D}i{\theta }_i{x}^{i-1}\right)}^2\mathrm{d}x$ with

$$m^{\prime }\left(x\right)=-\sum _{i=1}^k{w}_i\frac{x-{\mu }_i}{{\sigma }_i^2}p(x_i;{\mu }_i,{\sigma }_i),$$

denoting the derivative of the Gaussian mixture density $m\left(x\right)$. It follows that:

$$D_{H,2}[m:q]=\int m^{\prime }\left(x\right)\mathrm{d}x-2\sum _{i=1}^{D}i{\theta }_i\int x^{i-1}{m}^{\prime }\left(x\right)m\left(x\right)\mathrm{d}x+\sum _{i,j=1}^{D}ij{\theta }_i{\theta }_j\int x^{i+j-2}m{\left(x\right)}^2\mathrm{d}x,$$

where

$$\int x^i{m}^{\prime }\left(x\right)m\left(x\right)\mathrm{d}x=-\sum w_a{w}_b\int \frac{x-{\mu }_a}{{\sigma }_a^2}{x}^{i}p(x;{\mu }_a,{\sigma }_a)p(x;{\mu }_b,{\sigma }_b)\mathrm{d}x.$$

Therefore, we have

$$D_{H,2}[m:q]=\int m^{\prime }\left(x\right)\mathrm{d}x-2\sum _{i=1}^{D}i{\theta }_i\int x^{i-1}{m}^{\prime }\left(x\right)m\left(x\right)\mathrm{d}x+\sum _{i,j=1}^{D}ij{\theta }_i{\theta }_j\int x^{i+j-2}m{\left(x\right)}^2\mathrm{d}x$$

with $m^{\prime }\left(x\right)=-\sum w_a\frac{x-{\mu }_a}{{\sigma }_a^2}p(x;{\mu }_a,{\sigma }_a)$.

Since $p_a\left(x\right)p_b\left(x\right)={\kappa }_{a,b}p(x;{\mu }_{ab},{\sigma }_{ab})$, with

$$\begin{array}{ccc}\hfill {\mu }_{ab}& =& {\sigma }_a^2{\sigma }_b^2({\sigma }_b^2{\mu }_a+{\sigma }_a^2{\mu }_b),\hfill \\ \hfill {\sigma }_{ab}& =& \frac{{\sigma }_a{\sigma }_b}{\sqrt{{\sigma }_a^2+{\sigma }_b^2}},\hfill \\ \hfill {\kappa }_{a,b}& =& exp(F({\mu }_{ab},{\sigma }_{ab})-F({\mu }_a,{\sigma }_a)-F({\mu }_b,{\sigma }_b)),\hfill \end{array}$$

and

$$F(\mu ,\sigma )=\frac{{\mu }^2}{2{\sigma }^2}+\frac{1}{2}log\left(2\pi {\sigma }^2\right),$$

the log-normalizer of the Gaussian exponential family [42].

Therefore, we obtain

$$\int p_a\left(x\right)p_b\left(x\right)x^l\mathrm{d}x={\kappa }_{a,b}{m}_l({\mu }_{ab},{\sigma }_{ab}).$$

Thus, the Hyvärinen divergence $D_{H,2}$ of order 2 between a GMM and a PED is available in closed-form. □

For example, when $k=1$ (i.e., mixture m is a single Gaussian $p_{{\mu }_1,{\sigma }_1}$) and $p_{\theta }$ is a normal distribution (i.e., PED with $D=2$, $q_{\theta }=p_{{\mu }_2,{\sigma }_2}$), we obtain the following formula for the order-2 Hyvärinen divergence:

$$D_{H,2}[p_{{\mu }_1,{\sigma }_1}:p_{{\mu }_2,{\sigma }_2}]=\frac{{({\sigma }_1^2-{\sigma }_2^2)}^2+2{({\mu }_2-{\mu }_1)}^2{\sigma }_1^2}{8\sqrt{\pi }{\sigma }_1^3{\sigma }_2^4}.$$

4. Experiments: Jeffreys Divergence between Mixtures

In this section, we evaluate our heuristic to approximate the Jeffreys divergence between two mixtures m and $m^{\prime }$:

$${\tilde{D}}_J[m,m^{\prime }]:={({\overline{\theta }}_{\mathrm{SME}}\left(m^{\prime }\right)-{\overline{\theta }}_{\mathrm{SME}}\left(m\right))}^{\top }({\overline{\eta }}_{\mathrm{MLE}}\left(m^{\prime }\right)-{\overline{\eta }}_{\mathrm{MLE}}\left(m\right)).$$

Recall that stochastically estimating the JD between k-GMMs with Monte Carlo sampling using s samples (i.e., ${\widehat{D}}_{J,s}[m:m^{\prime }]$) requires $\tilde{O}\left(ks\right)$ and is not deterministic. That is, different MC runs yield fluctuating values that may be fairly different. In comparison, approximating $D_J$ by ${\tilde{D}}_J$ using ${\mathrm{\Delta }}_J$ by converting mixtures to D-order PEDs require $O\left(k{D}^2\right)$ time to compute the raw moments and $O\left(D^2\right)$ time to invert a Hankel moment matrix. Thus, by choosing $D=2k$, we obtain a deterministic $O\left(k^3\right)$ algorithm that is faster than the MC sampling when $k^2\ll s$. Since there are, at most, k modes for a k-GMM, we choose order $D=2k$ for the PEDs.

To obtain quantitative results on the performance of our heuristic ${\tilde{D}}_J$, we build random GMMs with k components as follows: $m\left(x\right)={\sum }_{i=1}^k{w}_i{p}_{{\mu }_i,{\sigma }_i}\left(x\right)$, where $w_i\sim U_i$, ${\mu }_i\sim -10+10U_1^{\prime }$ and ${\sigma }_i\sim 1+U_2^{\prime }$, where the $U_i$’s, and $U_1^{\prime }$ and $U_2^{\prime }$ are independent uniform distributions on $[0,1)$. The mixture weights are then normalized to sum up to one. For each value of k, we make 1000 trial experiments to gather statistics and use $s={10}^5$ for evaluating the Jeffreys divergence ${\widehat{D}}_J$ by Monte Carlo samplings. We denote by $\mathrm{error}:=\frac{|{\widehat{D}}_J-{\mathrm{\Delta }}_J|}{{\widehat{D}}_J}$ the error of an experiment.

Table 1. Comparison of ${\tilde{\mathrm{\Delta }}}_J(m_1,m_2)$ with ${\widehat{D}}_J(m_1,m_2)$ for random GMMs.

k	D	Average Error	Maximum Error	Speed-Up
2	4	0.1180799978221536	0.9491425404132259	2008.2323536011806
3	6	0.12533811294546526	1.9420608151988419	1010.4917042114389
4	8	0.10198448868508087	5.290871019594698	474.5135294829539
5	10	0.06336388579897352	3.8096955246161848	246.38780782640987
6	12	0.07145257192133717	1.0125283726458822	141.39097909641052
7	14	0.10538875853178625	0.8661463142793943	88.62985036546912
8	16	0.4150905507007969	0.4150905507007969	58.72277575395611

presents the results of the experiments for $D=2k$: The table displays the average error, the maximum error (minimum error is very close to zero, of order ${10}^{-5}$), and the speed-up obtained by our heuristic ${\mathrm{\Delta }}_J$. Those experiments were carried out on a Dell Inspiron 7472 laptop (equipped with an Intel(R) Core(TM) i5-8250U CPU at 1.60 GHz).

Notice that the quality of the approximations of ${\tilde{D}}_J$ depend on the number of modes of the GMMs. However, calculating the number of modes is difficult [43,70], even for simple cases [71,72].

Figure 4 displays several experiments of converting mixtures to pairs of PEDs to obtain approximations of the Jeffreys divergence.

Figure 5 illustrates the use of the order-2 Hyvärinen divergence $D_{H,2}$ to perform model selection for choosing the order of a PED.

Finally, Figure 6 displays some limitations of the GMM to PED conversion when the GMMs have many modes. In that case, running the conversion ${\overline{\eta }}_{\mathrm{MLE}}$ to obtain ${\tilde{\theta }}_T\left({\overline{\eta }}_{\mathrm{MLE}}\right)$ and estimate the Jeffreys divergence by

$${\tilde{\mathrm{\Delta }}}_J^{\mathrm{MLE}}[m_1,m_2]={({\tilde{\theta }}_2^{\mathrm{MLE}}-{\tilde{\theta }}_1^{\mathrm{MLE}})}^{\top }({\overline{\eta }}_2^{\mathrm{MLE}}-{\overline{\eta }}_1^{\mathrm{MLE}}),$$

improves the results but requires more computation.

Next, we consider learning a PED by converting a GMM derived itself from a Kernel Density Estimator (KDE). We use the duration of the eruption for the Old Faithful geyser in Yellowstone National Park (Wyoming, USA): The dataset consists of 272 observations (https://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat) (access date: 25 October 2021) and is included in the R language package ‘stats’. Figure 7 displays the GMMs obtained from the KDEs of the Old Faithful geyser dataset when choosing for each component $\sigma =0.05$ (left) and $\sigma =0.1$. Observe that the data are bimodal once the spurious modes (i.e., small bumps) are removed, as studied in [32]. Barron and Sheu [32] modeled that dataset using a bimodal PED of order $D=4$, i.e., a quartic distribution. We model it with a PED of order $D=10$ using the integral-based Score Matching method. Figure 8 displays the unnormalized bimodal density $q_1$ (i.e., ${\tilde{p}}_1$) that we obtained using the integral-based Score Matching method (with $\mathcal{X}=(0,1)$).

5. Conclusions and Perspectives

Many applications [7,73,74,75] require computing the Jeffreys divergence (the arithmetic symmetrization of the Kullback–Leibler divergence) between Gaussian mixture models. Since the Jeffreys divergence between GMMs is provably not available in closed-form [13], one often ends up implementing a costly Monte Carlo stochastic approximation of the Jeffreys divergence. In this paper, we first noticed the simple expression of the Jeffreys divergence between densities $p_{\theta }$ and $p_{{\theta }^{\prime }}$ of an exponential family using their dual natural and moment parameterizations [22] $p_{\theta }=p^{\eta }$ and $p_{{\theta }^{\prime }}=p^{{\eta }^{\prime }}$:

$$D_J[p_{\theta },p_{{\theta }^{\prime }}]={({\theta }^{\prime }-\theta )}^{\top }({\eta }^{\prime }-\eta ),$$

where $\eta =\nabla F\left(\theta \right)$ and ${\eta }^{\prime }=\nabla F\left({\theta }^{\prime }\right)$ for the cumulant function $F\left(\theta \right)$ of the exponential family. This led us to propose a simple and fast heuristic to approximate the Jeffreys divergence between Gaussian mixture models: First, convert a mixture m to a pair $(p_{{\overline{\theta }}^{\mathrm{SME}}},p^{{\overline{\eta }}^{\mathrm{MLE}}})$ of dually parameterized polynomial exponential densities using extensions of the Maximum Likelihood and Score Matching Estimators (Theorems 1 and 3), and then approximate the JD deterministically by

$$D_J[m_1,m_2]\simeq {\tilde{D}}_J[m_1,m_2]={({\tilde{\theta }}_2^{\mathrm{MLE}}-{\tilde{\theta }}_1^{\mathrm{MLE}})}^{\top }({\overline{\eta }}_2^{\mathrm{MLE}}-{\overline{\eta }}_1^{\mathrm{MLE}}).$$

The order of the polynomial exponential family may be either prescribed or selected using the order-2 Hyvärinen divergence, which evaluates in closed form the dissimilarity between a GMM and a density of an exponential-polynomial family (Theorem 4). We experimentally demonstrated that the Jeffreys divergence between GMMs can be reasonably well approximated by ${\tilde{D}}_J$ for mixtures with a small number of modes, and we obtained an overall speed-up of several order of magnitudes compared to the Monte Carlo sampling method. We also propose another deterministic heuristic to estimate $D_J$ as

$${\tilde{D}}_J^{\mathrm{MLE}}[m_1:m_2]={({\tilde{\theta }}_2^{\mathrm{MLE}}-{\tilde{\theta }}_1^{\mathrm{MLE}})}^{\top }({\overline{\eta }}_2^{\mathrm{MLE}}-{\overline{\eta }}_1^{\mathrm{MLE}}),$$

where ${\tilde{\theta }}^{\mathrm{MLE}}\approx \nabla F\left({\overline{\eta }}^{\mathrm{MLE}}\right)$ is numerically calculated using an iterative conversion procedure based on maximum entropy [58] (Section 2.3.1). Our technique extends to other univariate mixtures of exponential families (e.g., mixtures of Rayleigh distributions, mixtures of Gamma distributions, or mixtures of Beta distributions, etc). One limitation of our method is that the PED modeling of a GMM may not guarantee obtaining the same number of modes as the GMM even when we increase the order D of the exponential-polynomial densities. This case is illustrated in Figure 9 (right).

Although PEDs are well-suited to calculate Jeffreys divergence compared to GMMs, we point out that GMMs are better suited for sampling, while PEDs require Monte Carlo methods (e.g., adaptive rejection sampling or MCMC methods [62]). Furthermore, we can estimate the Kullback–Leibler divergence between two PEDs using rejection sampling (or other McMC methods [62]) or by using the $\gamma$-divergence [76] with $\gamma$ close to zero [66] (e.g., $\gamma =0.001$). The web page of the project is https://franknielsen.github.io/JeffreysDivergenceGMMPEF/index.html (accessed on 25 October 2021).

This work opens up several perspectives for future research: For example, we may consider bivariate polynomial-exponential densities for modeling bivariate Gaussian mixture models [29], or we may consider truncating the GMMs in order to avoid tail phenomena when converting GMMs to PEDs [77,78].