An Efficient Gaussian Mixture Model and Its Application to Neural Network

1. Introduction

Gaussian Mixture Models (GMMs) are powerful, parametric, and flexible models in probabilistic modeling. In particular, GMMs address two important problems: the inverse problem and the uncertainity problem. The first, known as the inverse problem, arises when $f^{-1}\left(x\right)$ and x are not unique one-to-one mappings [23]. The second issue involves the uncertainty of both the data (aleatoric) and the model (epistemic) [4,5,6,7]. When considering these two problems together, the classic modeling assumption $y=f\left(x\right)+ϵ$ becomes inadequate. In a probabilistic sense, this assumption is equivalent to $y\sim \mathit{Normal}\left(f\left(x\right),\sigma \right)$. However, to address the inverse and uncertainty problems, the model assumption should be $y\sim \mathit{UnkownDistribution}\left(\theta \left(x\right)\right)$. Under these conditions, when y no longer follows any known distribution, GMM become a viable option for modeling the data.

Fitting a GMM to observed data is related to two problems: determining the mixture components needed, and estimating the parameters of each mixture component. For the first problem, various techniques have been introduced [1]. For the second problem, the expectation maximization (EM) algorithm [11], developed as early as 1977, can fit GMMs by minimizing the negative log-likelihood (NLL). While the EM algorithm remains one of the most popular choices for learning GMMs [10,12,13], other methods have also been developed. However, the EM algorithm has some drawbacks. Srebro [16] raised an open question regarding the existence of local maxima between an arbitrary distribution and a GMM. This study also addressed the possibility that the Kullback-Leibler (KL) divergence between an arbitrary distribution and mixture models may have non-global local minima, especially when the target distribution has more than k components, but we fit it with a GMM with fewer than k components. Améndola et al. [15] demonstrated that the likelihood function for GMMs is transcendental. Jin et al. [14] resolved the problem Srebro discovered and proved that with random initialization, the EM algorithm will converge to a bad critical point with high probability, and the local maxima can be arbitrarily worse than any global optimum. The maximum likelihood function is transcendental, resulting in multiple critical points. Because of the problems related to likelihood functions, some methods, such as the Wasserstein distance, have been explored. Wasserstein distances are presented to replace the NLL function, which minimizes either the KL divergence or the NLL function [10,17,18]. These methods avoid the disadvantage of using NLL; however, the calculation and formulation of the learning process are relatively complex. Another state-of-the-art method is Bayesian variational inference (BVI) for a Gaussian mixture model [33,35,37]. This method applies Bayesian techniques and variational inference to the Gaussian mixture model. Practically, this method avoids the singularity problems in EM; however, many priors are required in the inference algorithms. When these priors are mismatched, the algorithm’s performance may be less accurate than EM.

Instead of following the EM and BVI, we adopt a different path. One of the main contributions of this study is the concept of GMM expansion, which is based on the idea of Fourier expansion applied to density approximations. Like Fourier expansion, we process that any density can be approximated by a linear combination of a set of Gaussian distributions. If this argument is valid, the Gaussian distributions in the GMM become the basis, similar to the cosine and sine functions in Fourier expansions. Some studies show that all normal mixture densities are dense in the set of all density functions under the $L^1$ metric [1,2,3,34], and the mixture densities are dense in the set of all density functions under the Wasserstein metric. In this study, we further prove that mixture densities built under the GMM expansion can approximate an arbitrary density with a certain accuracy. Numerical analysis was also conducted to secure the usability of the GMM expansion in practice.

In GMM expansion, the Gaussian distributions act as "basis functions" that do not require additional parameterization. The only parameter that needs to be learned is $\pi$. These "basis functions" represent how much detail our models are able to capture. The benefit of this idea is that learning on a predetermined basis becomes relatively simple. A simple algorithm without heavy formulation is proposed for density approximation under the GMM expansion. We show that our method can be more accurate, faster, and easier to use compared to state-of-the-art methods such as EMs and BVI.

GMMs have been used in computer vision and pattern classification, among many other fields [10,19,20,21,22]. Another significant benefit of the proposed method is that the GMM expansion can improve the integration of GMM in neural networks. The proposed method has no singularity issue compared to EM, which significantly improves performance when a neural network is trained to predict a distribution. The simplicity of the proposed method also allows us to easily utilize GMM into neural networks to manipulate the latent space. We built three applications to demonstrate the benefits that GMM can bring to neural networks.

Section 2 introduces the GMM expansion method, which proves that GMM expansion can approximates arbitrary densities, learning algorithms, and convergence. Section 3 provides our numerical experiments and a comparison with state-of-the-art density learning methods. Section 4 presents three neural-network applications that demonstrate the integration of the proposed method into a neural network. The final section provides a summary and conclusions.

2. GMM Expansion

In this section, we provide a detailed description of the principle of the GMM expansion using learning algorithms. Theoretical and numerical analyses were conducted to prove that the GMM expansion can approximate arbitrary densities with a bounded error. Gaussian mixture models in the high-dimensional case are theoretically possible; however, similar to the Fourier expansion method, the complexity of the algorithm for GMM grows exponentially with the number of dimensions. Due to arithmetic constraints, we will discuss higher dimensions in future work. In this study, we conducted experiments in one and two dimensions.

2.1. The Concept of GMM Expansion

The idea behind the GMM decomposition is akin to a Fourier series in probability, where an arbitrary density can be approximately decomposed by an infinite mixture of Gaussian distributions. As aforementioned, several studies have proven that the likelihood of a GMM has numerous critical points. There are multiple sets of $\pi ,\mu$ and $\sigma$ that can produce the same likelihood result on a particular dataset. If we consider that GMM can decompose any density like the Fourier series decomposes into a periodic sequence, we can hold $\mu$ and $\sigma$ unchanged and focus on finding the optimal solution for $\pi$. Based on this concept, the setup of the GMM expansion is as follows:

• Observation X follows an unknown distribution with density $f\left(x\right)$;
• $g\left(x\right)={\sum }_N{\pi }_i{\varphi }_i\left(x\right)$ is the density of GMM;
• ${\sum }_N{\pi }_i=1$;
• ${\varphi }_i\left(x\right)$ is the i-th normal distribution with mean ${\mu }_i$ and standard deviation $\sigma$;
• define $r={\displaystyle \frac{\max \left(X\right)-\min \left(X\right)}{n}}$ as the interval for locating ${\mu }_i$;
• ${\mu }_i=\min \left(X\right)+i\times r$ ;
• $\sigma$ is a hyperparameter and the same for all Gaussian components, usually $r\le \sigma \le 5r$.

The $\mu$ and $\sigma$ of the component distributions are non-parametric and do not need to be optimized. The parameter $\pi$ is the only remaining problem for us to solve. There are justification $\sigma$ as a hyper-parameter, which we will further demonstrate in the next subsection. From the above settings, to ensure this system works, the GMM expansion must satisfy one condition: for any $f\left(x\right)$, there exists a $g\left(x\right)$ such that $\mid f\left(x\right)-g\left(x\right)\mid <ϵ$. In the next subsection, we will prove that this is indeed the case.

2.2. Approximate Arbitrary Density by GMM Expansion

First, we would like to introduce a Lemma from Villani’s book [34].

Lemma 1 

[34] Let $\mathcal{X}$ be a complete separable metric space and $P\in (0,\infty )$. Then the Wasserstein space $P_p\left(\mathcal{X}\right)$, metrized by the Wasserstein distance $W_p$, is also a complete separable metric space. i.e. any probability measure can be approximated by a sequence of probability measures with finite support.

The proof of this Lemma is given by Theorem 6.18 in [34]. Under this Lemma, the Dirac measure, which can be seen as a degenerate GMM, is able to approximate any arbitrary distribution with arbitrary precision in the Wasserstein space. To further elaborate on this, when the Gaussians are very close to degeneration, which is close to the Dirac measure, we have the following numerical interpretation of Lemma 1 under the GMM expansion setting.

Consider a fine enough partition based on $\mu$ over the support of the probability measure such that ${\omega }_i=\left({\mu }_i-r,{\mu }_i+r\right]$, then we have:

$$\left[\begin{array}{cccc}{\int }_{{\omega }_0}{\varphi }_0\left(x\right)dx& .& .& .{\int }_{{\omega }_0}{\varphi }_n\left(x\right)dx\\ ⋮& ⋮& ⋮& ⋮\\ {\int }_{{\omega }_N}{\varphi }_0\left(x\right)dx& .& .& .{\int }_{{\omega }_N}{\varphi }_n\left(x\right)dx\end{array}\right]\left[\begin{array}{c}{\pi }_0\\ ⋮\\ {\pi }_N\end{array}\right]=\left[\begin{array}{c}{P}_g\left\{{\omega }_0\right\}\\ ⋮\\ P_g\left\{{\omega }_N\right\}\end{array}\right].$$

If ${\int }_{{\omega }_i}{\varphi }_i\left(x\right)dx\approx 1$, indicating that each Gaussian component is very dense in its region, the GMM becomes a discrete distribution that can precisely map the estimated probability mass. It becomes:

$$\left[\begin{array}{cccc}1& .& .& .0\\ ⋮& ⋮& ⋮& ⋮\\ 0& .& .& .1\end{array}\right]\left[\begin{array}{c}{\pi }_0\\ ⋮\\ {\pi }_N\end{array}\right]\approx \left[\begin{array}{c}{P}_g\left\{{\omega }_0\right\}\\ ⋮\\ P_g\left\{{\omega }_N\right\}\end{array}\right]$$

$$\left[\begin{array}{c}{P}_g\left\{{\omega }_0\right\}\\ ⋮\\ P_g\left\{{\omega }_N\right\}\end{array}\right]\approx \left[\begin{array}{c}{\pi }_0\\ ⋮\\ {\pi }_N\end{array}\right]=\left[\begin{array}{c}{P}_f\left\{{\omega }_0\right\}\\ ⋮\\ P_f\left\{{\omega }_N\right\}\end{array}\right]$$

(1)

where $P_g$ is the density function for the GMM and $P_f$ is the target density function. These simple deductions, together with Lemma 1 show that any probability measure can also be approximated by probability measures constructed by the GMM.

Theorem 2. 

For any probability density $f\left(x\right)$ and a given error e, we can find a Gaussian mixture $g\left(x\right)$ to approximate it. Proof Let $f\left(x\right)$ be a piecewise continuous density function, whose integral over a finite interval $[a,b]$ is

$$\begin{array}{c}\hfill {\int }_a^{b}f\left(x\right)dx=1-ϵ,ϵ>0.\end{array}$$

(2)

We consider a sequence $a=x_1<x_2<\cdots <x_{N-1}<x_N=b$ such that $f\left(x\right)$ is continuous in every interval $[x_i,x_{i+1}]$ for $x_i\le x<x_{i+1}$ and $i=1,2,\cdots ,N-1$. When the partitioning is sufficiently fine, continuous function $f\left(x\right)$ on a bounded interval $[x_i,x_{i+1}]$ is bounded and there exist a constant ${\xi }_i$ such that

$$\begin{array}{c}\hfill |f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|<{\xi }_i.\end{array}$$

(3)

The error of two functions under one norm is

$$\begin{array}{c}\hfill \sum _{i=1}^N|f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|<\sum _{i=1}^N{\xi }_i,\end{array}$$

(4)

and

$$\begin{array}{cc}& |{\sum }_{i=1}^{N}f\left(x_i\right)(x_{i+1}-x_i)-{\int }_a^{b}f\left(x\right)dx|\end{array}$$

(5)

$$\begin{array}{cc}& \le {\sum }_{i=1}^N|f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|.\end{array}$$

(6)

There exists a sufficiently large $\widehat{N}$ such that when $N>\widehat{N}$, for any $\xi$

$$\begin{array}{c}\hfill \sum _{i=1}^N|f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|<\sum _{i=1}^N{\xi }_i<\xi .\end{array}$$

(7)

According to Lemma 1, the first term of Equation (6) can be approximated with an arbitrary precision using a finite combination of Dirac measures. We assume that the error between the Gaussian components ${\Phi }_i\left(x\right)$ and the Dirac measure $\delta \left(x\right)$ within each interval is $e_i$ in one norm.

$$\begin{array}{c}\hfill f\left(x_i\right)(x_{i+1}-x_i)\approx {\alpha }_i\delta (x;[x_i,x_{i+1}])\approx {\alpha }_i{\Phi }_i\left(x\right)\end{array}$$

(8)

where $0\le {\alpha }_i\le 1,\sum {\alpha }_i\approx 1$.

The global error can be shown as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \sum _{i=1}^N|{\alpha }_i\Phi (x;[x_i,x_{i+1}],\sigma )-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|\\ \hfill =\sum _{i=1}^N|{\alpha }_i\Phi (x;[x_i,x_{i+1}],\sigma )-f\left(x_i\right)(x_{i+1}-x_i)\\ \hfill +f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|\\ \hfill \le \sum _{i=1}^N|{\alpha }_i{\Phi }_i(x;[x_i,x_{i+1}],\sigma )-f\left(x_i\right)(x_{i+1}-x_i)|\\ \hfill +\sum _{i=1}^N|f\left(x_i\right)(x_{i+1}-x_i)-{\int }_{x_i}^{x_{i+1}}f\left(x\right)dx|\\ \hfill <\sum _{i=1}^N{e}_i+\sum _{i=1}^N{\xi }_i\\ \hfill <\sum _{i=1}^N{e}_i+\xi \end{array}\end{array}$$

(9)

and

$$\begin{array}{c}\hfill |\sum _{i=1}^N{\alpha }_i{\Phi }_i(x;[a,b],\sigma )-{\int }_{-\infty }^{\infty }f\left(x\right)dx|<\sum _{i=1}^N{e}_i+\xi +ϵ=e.\end{array}$$

(10)

where ${\sum }_{i=1}^N{\alpha }_i{\Phi }_i(x;[a,b],\sigma )$ is a mixture of Gaussian.

Remark 1 

Despite the introduction of multiple error terms in the above approximation, the global error can be controlled, and the results of the numerical experiments support this conclusion. The GMM expansion, as a special kind of Gaussian mixture model, also satisfies the above theorem, and an example is given in the Appendix A.

2.3. Total Variation Distance Analysis For Small Region

Even though GMMs can approximate arbitrary densities similar to frequency distributions under Wasserstein metrics, the density within a fine region remains undiscussed. In other words, even if we learn a model where the probability mass is equal to the target density $P_g\left\{\omega \right\}=P_f\left\{\omega \right\}$, the differentiation might not be equal ${\displaystyle \frac{\mathrm{d}{P}_g\left\{\omega \right\}}{\mathrm{d}x}}\ne {\displaystyle \frac{\mathrm{d}{P}_f\left\{\omega \right\}}{\mathrm{d}x}}$. To further analyze the difference in probability density in fine regions, we introduce the well-known total variation distance.

Assuming we have learned a Gaussian mixture density g that correctly maps to the target distribution fand considering the same partition as shown in Section 2.2, let Equation (1) be satisfied. The total variation distance (TVD) based on the partition $\omega$ is given by

$$TVD={\displaystyle \frac{1}{2}}\sum _{i=1}^n{\int }_{{\omega }_i}|f\left(x\right)-g\left(x\right)|dx.$$

Because the integration is difficult to calculate, we divide ${\omega }_i$ further into m smaller regions to estimate TVD:

$$\sum _{k=1}^m{d}_k={\omega }_i,$$

$$\begin{array}{cc}\hfill TVD& \approx {\displaystyle \frac{1}{2}}\sum _{i=1}^n\sum _{j=1}^m{\int }_{d_{ij}}|f\left(x\right)-g\left(x\right)|dx\hfill \\ & \approx {\displaystyle \frac{1}{2}}\sum _{i=1}^n\sum _{j=1}^m|P_f\left\{d_{ij}\right\}-P_g\left\{d_{ij}\right\}|.\hfill \end{array}$$

Denote that

$$L\left({\omega }_i\right)={\displaystyle \frac{1}{2}}\sum _{j=1}^m|P_f\left\{d_{ij}\right\}-P_g\left\{d_{ij}\right\}|.$$

(11)

Here $L\left({\omega }_i\right)$ is the discrete approximation of TVD within region ${\omega }_i$. If $L\left({\omega }_i\right)$ is universally bounded in some conditions for all $\omega$, the overall TVD is bounded. For any $x_{ij}$ from ${\omega }_i$, if

$$|P_f\left\{x_{ij}\pm d\right\}-P_g\left\{x_{ij}\pm d\right\}|\le {\tau }_i\le \tau ,$$

we obtain

$$TVD\le n\ast m\ast \tau .$$

Remark 2 In the Appendix B, we illustrate why TVD is controllable and $L\left({\omega }_i\right)$ is universally bounded under certain conditions. The experiments described in Section 3 verify this result.

2.4. Learning Algorithm and Convergence

From Section 2.1 to Section 2.3, we draw the following conclusion: The parameter $\pi$s are proportional to the distribution of data, resulting in a simple update for the learning algorithm. The whole algorithm can be demonstrated as follows.

Algorithm 1:GMM Expansion

Initialization Calculate r based on the number of Gaussian components n, $r={\displaystyle \frac{\mathrm{Max}\left(X\right)-\mathrm{Min}\left(X\right)}{n}}$. Evenly spread $\mu =\left[{\mu }_1,{\mu }_2,\cdots ,{\mu }_n\right]$ across $\mathrm{Max}\left(X\right)$ and $\mathrm{Min}\left(X\right)$, so that ${\mu }_{i+1}-{\mu }_i=r$. Define hyperparameter $\sigma$ with respect to r, e.g., $\sigma =r$, $\sigma =3r$, etc. Define hyperparameter d with respect to $\sigma$, e.g., $d=\sigma /4$, $d=\sigma /6$, etc. Define a constant $\alpha =\beta (P\left({\mu }_i\pm {\displaystyle \frac{d}{2}}\right)-P\left({\mu }_i-{\displaystyle \frac{r}{2}}\le x\le {\mu }_i-{\displaystyle \frac{r}{2}}+d\right))$ . Initialize ${\pi }_1={\pi }_2=\cdots ={\pi }_N={\displaystyle \frac{1}{N}}$. Update Procedure 1: for each $x_d$ do 2:     $j=\mathrm{argmin}(\mid x_d-\mu \mid )$ 3:     ${\pi }_j^{+1}={\pi }_j+\alpha$ 4: end for Finally, re-normalize $\pi$s to ensure $\sum {\pi }_i=1$.

Since we always seek to maximize the likelihood, a gradient ascent method should be applied. The derivative of the i-th $\pi$ with respect to$g\left(x_d\right)$ for each observation partially reveals the methodology of the proposed algorithm:

$${\displaystyle \frac{\partial g\left(x_d\right)}{\partial {\pi }_i}}={\varphi }_i\left(x_d\right).$$

(12)

Because $\mu$, and $\sigma$ are not longer parameters, the gradient direction from every single data point is always positive, which is directly related to $\varphi$. After learning from the whole dataset, the final step is re-normalization to ensure $\sum {\pi }_i=1$. Figure 2 illustrates the learning process. The parameter $\sigma$ is regarded as a hyperparameter that controls the smoothness of the GMM. More details are provided in Appendix B. The parameter $\alpha$ is a proportional hyper-parameter designed based on the proof shown in Section 2.1, Section 2.2 and Section 2.3. The parameter $\beta$ is the learning rate, similar to any gradient ascend method. The probability part of $\alpha$ can be referred to in Eq. A7, where P is the probability of a normal distribution. The reasoning is provided in the Appendix B and illustrated in Figure 2. Considering a data point, the algorithm adjusts the probability with respect to $\pi$ while maintaining the rest of the density in proportion. Further explanation is provided in the convergence analysis. The proposed algorithm is simple yet efficient; it processes each data point in the dataset only once. On the contrary, EM and BVI algorithms require several iterations over the dataset to converge.

Moving on to the discussion of the convergence, consider that we have a set of bases Gaussian components ${\left\{{\varphi }_i\left(x\right)\right\}}_n$. According to Theorem 6.18, the mixture densities are dense and can approximate any Dirac mass with arbitrary precision under the Wasserstein distance. Dirac masses can be considered a degeneration of the Gaussian distribution. Assuming the Gaussian components are sufficiently dense, for any sample data point $x_d$, we can find a Gaussian component very close to it, such that $\varphi (x_d;{\mu }_d,\sigma )\approx \varphi ({\mu }_d;{\mu }_d,\sigma )$. This means:

$$\begin{array}{c}g\left(x_d\right)=\sum _i^N{\pi }_i{\varphi }_i\left(x_d\right)\approx {\pi }_d{\varphi }_d\left(x_d\right)\approx {\pi }_d{\varphi }_d\left({\mu }_d\right),\end{array}$$

(13)

$$\begin{array}{c}{\varphi }_d\left({\mu }_d\right)=a,g\left(x_d\right)\approx a{\pi }_d,\end{array}$$

(14)

where a is a constant. Similar to BVI and the EM algorithm, the proof of convergence and the algorithm both rely on a binary latent variable Z[33,35,36]. In this context, we assume $\gamma$ is a binary latent variable, which indicates that observation $x_d$ is sampled from the n-th mixture component. This ensures that the density function is expected as:

$$\begin{array}{c}p\left(x\right|\gamma )=\prod _{n=1}^N\varphi {\left(x\right|{\mu }_n,\sigma )}^{\gamma },\end{array}$$

(15)

We can observe the connection between the assumptions we based on the analyses provided in [34] and the latent variable model used in BVI and EM. This indicates that Eq.(13)is equivalent to Eq.(15). Thus, we will pursue a similar approach by considering a variable $Z\in \mathbb{N}={\left\{z_n\right\}}_N$.

$$\begin{array}{c}{z}_i=\sum _{d=1}^D{\gamma }_{di},\\ \sum _{i=1}^N{z}_i=D,\\ {\gamma }_{di}=\left\{\begin{array}{cc}1\hfill & i=\mathrm{argmin}(x_d-\mu )\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

The log-likelihood function $log\left(\mathcal{L}\right)$ is:

$$\begin{array}{c}log\left(\mathcal{L}\right)=\sum _d^{D}log\left(\sum _i^N{\pi }_i{\varphi }_i\left(x_d\right)\right),\\ =\sum _d^{D}log\left(\sum _i^N{\gamma }_{di}{\pi }_i{\varphi }_i\left(x_d\right)\right),\\ =\sum _d^D\sum _i^N{\gamma }_{di}log\left({\pi }_i{\varphi }_i\left(x_d\right)\right),\\ =\sum _i^N\sum _d^D{\gamma }_{di}log\left({\pi }_i{\varphi }_i\left(x_d\right)\right),\\ log\left(\mathcal{L}\right)\approx \sum _i^N{z}_{i}log\left(a{\pi }_i\right).\end{array}$$

Define a function with the Lagrange multiplier :

$$\begin{array}{c}\mathcal{A}=\sum _i^N{z}_{i}log\left(a{\pi }_i\right)+\lambda (1-\sum _i^N{\pi }_i)\end{array}$$

$$\begin{array}{c}{\displaystyle \frac{\partial \mathcal{A}}{\partial {\pi }_i}}=z_i{\displaystyle \frac{1}{{\pi }_i}}-\lambda =0,\\ \mathrm{Set}\lambda =D,\mathrm{then}{\pi }_i^{*}={\displaystyle \frac{z_i}{D}}\end{array}$$

From Algorithm 1, denote ${\pi }^{alg}$ is the $\pi$ learned from the proposed algorithm. After running the update procedure, we obtain:

$$\begin{array}{c}{\pi }_i^{alg}={\displaystyle \frac{\frac{1}{N}+z_i\alpha }{{\sum }_i^N\frac{1}{N}+z_i\alpha }}={\displaystyle \frac{\frac{1}{N}+z_i\alpha }{1+D\alpha }},\end{array}$$

when N and $D\alpha$ are sufficiently large,

$$\begin{array}{c}{\pi }_i^{alg}\approx {\pi }^{*},\end{array}$$

Based on the above assumption, our method is equivalent to the statistic of the frequency distribution.

$$\begin{array}{c}g\left(x_d\right|\pi )=\sum _i{\displaystyle \frac{1}{N}}\varphi \left(x_d\right)\approx {\displaystyle \frac{1}{N}}{\varphi }_k\left(x_d\right)\\ g\left(x_d\right|{\pi }^{*})=\sum _i{\displaystyle \frac{z_i}{D}}\varphi \left(x_d\right)\ge {\displaystyle \frac{z_k}{D}}{\varphi }_k\left(x_d\right)\\ z_i\ge 1,D\le N,k=\mathrm{argmin}(x_d-\mu ),\end{array}$$

For every data point $x_d$,

$$\begin{array}{c}g\left(x_d\right|{\pi }^{*})>g\left(x_d\right|\pi ),\\ \mathrm{Likelihood}\mathcal{L}\left(X\right|{\pi }^{alg})\approx \mathcal{L}\left(X\right|{\pi }^{*})\ge \mathcal{L}\left(X\right|\pi )\end{array}$$

The proposed algorithm can employ both ${\pi }^{alg}$ and ${\pi }^{*}$ for approximation. ${\pi }^{alg}$ is ideal for situations with a small number of data points where a conservative approximation for unseen regions is required. On the other hand, ${\pi }^{*}$ is more suitable when a closer approximation to the dataset is needed.

3. Numerical Experiments

In this section, we demonstrate the effectiveness of the proposed GMM expansion and learning algorithm. Section 3.1 presents experiments that illustrate and discuss the fundamentals of the proposed method for approximating arbitrary densities. For usability, Section 3.2 provides comparative studies, contrasting the proposed method with state-of-the-art techniques such as BVI and EM algorithm. These comparisons highlight the advantages of the proposed method.

3.1. Density Approximation

The GMM expansion aims to approximate densities without relying on the specific form of the target distributions. Figure 1 shows cases of two density approximation. The left side of Figure 1 displays the target densities, while the right side shows these target densities alongside their approximations using our GMM method. In the upper plot, the target density is a mixture of continuous and smooth distributions, using five normal distributions. The lower plot features a target distribution that combines normal, T, and uniform distributions.

In terms of approximation accuracy, the proposed method performed exceptionally well for both the smooth and non-smooth cases. The overall and most crucial features of the target densities were well-preserved, although some details of non-smooth densities were difficult to recover. Considering that the Gaussian distributions are inherently smooth, we expected the GMM to perform better when approximating smooth target densities. Geometrically, when target distributions consist of components with bell-shaped structures, GMM is a particularly suitable choice for approximation. However, in cases like the target distribution in the lower plot, which includes a uniform distribution among its components, GMM may struggle to capture fine details unless a sufficiently large number of components are used. According to the GMM expansion idea, GMM is expected to capture finer structures, such as deep cuts and straight lines, with a larger set of basis components. This expectation is consistent with the theoretical proof provided in Section 2. As long as GMM captures most of the critical features with a certain level of accuracy, the approximation can be considered good. Further numerical justification is provided later in this section. Further numerical justification is provided later in this section.

The learning process is illustrated in Figure 2. Initially, the density of our GMM closely resembles a uniform distribution before training. As the training progresses and the dataset is processed, the density gradually conforms to the original target density. Given that the base distributions are fine and dense, we can achieve a good approximation without needing to learn the parameters $\mu$ and $\sigma$ for each Gaussian distribution. Unlike traditional methods where GMM requires learning $\pi$$\mu$ and $\sigma$, fixing $\mu$ and $\sigma$ as bases and smoothness, especially in neural network, which we will present in section 4 adjusting the density through $\pi$ offers advantages, particularly in neural network applications, as discussed in Section 4. Additional approximation examples are shown in Appendix Figure A5.

After graphically demonstrating the GMM expansion and learning algorithm, a quantitative study was conducted to further discuss the properties of the proposed method. In these experiments, the total variance distance(TVD) from Section 2 was applied. Since $\mu$ was predetermined, TVD is denoted as L and is calculated as follows.

$$\begin{array}{cc}\hfill 2L=& |P_{target}\left\{x\le {\mu }_1-d\right\}-P_{gmm}\left\{x\le {\mu }_1-d\right\}|\hfill \\ \hfill & +\sum _{i=1}^n|P_{target}\left\{{\mu }_i\pm d\right\}-P_{gmm}\left\{{\mu }_i\pm d\right\}|\hfill \\ \hfill & +|P_{target}\left\{x>{\mu }_N+d\right\}-P_{gmm}\left\{x>{\mu }_N+d\right\}|.\hfill \end{array}$$

(16)

Measuring the estimation accuracy using Equation (16) has several advantages. First, it is easy to calculate. Second, even without the target density in its current form, we can still compute the loss using a discrete statistical estimation. The minimum and maximum values of the loss, which were $[0,2]$, are clearly evident. Evidence is provided below.

$$\begin{array}{cc}\hfill 2L& \le |P_{target}\left\{x\le {\mu }_1-d\right\}-0|+\sum _{i=1}^n|P_{target}\left\{{\mu }_i\pm d\right\}-0|\hfill \\ \hfill & +|P_{target}\left\{x>{\mu }_N+d\right\}-0|+|0-P_{gmm}\left\{x\le {\mu }_1-d\right\}|\hfill \\ \hfill & +\sum _{i=1}^n|0-P_{gmm}\left\{{\mu }_i\pm d\right\}|+|0-P_{gmm}\left\{x>{\mu }_N+d\right\}|.\hfill \end{array}$$

That means $0\le L\le 1$.

In this quantitative study, the proposed method was used to approximate two sets of randomly generated distributions, calculating the average TVD using Equation (16) between the target distribution and the approximated GMM. Table 1 and Table 2 present the results of using our method to learn two types of distributions. Similar to Figure 1, the target distributions in Table 1 are randomly generated Gaussian mixture distributions. We sampled 5000 data points for learning and used GMM with 200 components for learning. More approximation examples are shown in Appendix Figure A5 subplots (1)–(3). Table 2 presents the results of approximating target distributions containing Normal, T, and uniform distributions, with examples shown in Appendix Figure A5 subplots (4)–(6). Since $0\le L\le 1$, the average L for both types of distributions are less than $0.1$. This indicates that approximation error is considerably low, with the margin of error for our experiments being less than $5\%$. In the later part of the section, we present a direct comparison with other methods. The results also indicate that the GMM generally performs better in approximating smooth target densities. Another important factor in the algorithm is the hyperparameter $\sigma$. Table 1 and Table 2 show how changes in $\sigma$ affect the accuracy of the algorithm. The approximation accuracy benefits from a slightly larger $\sigma$ for smooth densities, but a smaller $\sigma$ works better for non-smooth densities. However, the overall improvement in L function was less than $0.02$, indicating that the benefits are not significant. This result supports treating $\sigma$ as a hyperparameter. Depending on the situation, $\sigma$ can be set smaller (e.g., $1\ast r$) if a closer match to the original dataset is required or larger (e.g., $3\ast r$) for smoother distributions.

The above experiment discusses approximation accuracy with the same component distribution sizes but varying $\sigma$. It is important to address how the sizes of the component distributions and the dataset affect the approximation accuracy. The experimental results shown in Table 3 demonstrates the performance of the proposed method under various setups. Similar to the experiments in Table 2 and Table 1, the average L over approximations of 30 different randomly generated mixture densities was calculated. A GMM with component sizes ranging from 10 to 1000 was tested. Table 1, Table 2, and Table 3 show how increasing the data size from 5,000 to 50,000 affects approximation accuracy. The findings indicate that the proposed method benefits from a large dataset. Results from Table 1 and Table 2 show that with 5,000 data points and 200 components, the approximation accuracy was improved from 0.04369 to 0.01475 when the data size increased to 20,000. Further increasing the number of data points to 50,000, provided only a slight additional benefit. Table 3 indicates increasing the component size from 10 to 200 delivers significant improvements. However, increasing the component size to 500 and 1000 does not necessarily result in better performance. An approximation with 500 components and 50,000 data points achieved the best accuracy with $L=0.0101$, but the improvement over 200 components ($L=0.0125$) was minimal. These findings align with our intuitive understanding of GMM expansion, confirming that a certain number of basis components is necessary to ensure a good approximation. Increasing data points and components helps capture finer details, but determining the optimal number of components relative to the dataset size lacks a clear formula. Nonetheless, our experiments show that the proposed method is robust with 200 components.

The aforementioned experiments were conducted using one-dimensional settings. The proposed method as also applicable in high-dimensional cases. A two dimensional example is shown in Figure 3. In Fourier expansion, increasing dimensions exponentially enlarges the size of the base frequency components. The proposed method faces a similar problem, managing higher-dimensions requires more bases to capture inter-dimensional correlations. To address this, clustering methods such as K-means are used with the EM algorithm, where cluster sizes are predetermined. The BVI method attempts to solve this problem by optimizing the posterior distribution through techniques like mean-field approximation or Markov chain monte-carlo (MCMC). However, these methods each have limitations. For example, the EM algorithm can suffer from singularities, and the BVI approach requires prior and posterior distributions, which introduces implicit biases. Higher-dimensional density estimation is more complex compared to existing methods due to the correlation between dimensions. We plan to further develop this method and discuss future work. In the next subsection, we will directly compare the proposed method with state-of-the-art distribution approximation techniques, particularly EM and BVI.

3.2. Comparison Study

EM and BVI are state-of-the-art density approximation methods [33,35,37]. Scikit-learn1 provides a well-built Python package for implementing these methods, alongwith detailed explanation and summaries. In this subsection, we compare our proposed method with EM and BVI. The algorithms for both the EM and BVI were obtained from Scikit-learn.

Figure 4. Experiment Density Set

Figure 4. Experiment Density Set

The comprehensive performance comparison results for the proposed method, EM, and BVI are presented in Table 4. Accuracy and training time were compared, with all values presented as mean and standard deviation. Two variations of the EM algorithm were tested: EM K-means, which initializes parameters using the K-means methods, and EM-Random, which uses randomly selected data points. Similar to the experiments in Table 1, Table 2, and Table 3, a set of mixture distributions was randomly generated as illustrated in Figure 4. For each mixture distribution, 2000 data points were sampled for training. The component sizes ranged from 20 to 1000. We used the Total Variation Distance (TVD) to evaluate approximation accuracy. Since the cumulative density function is not provided in Scikit-learn, we estimated the TVD using a slightly modified approach. The equation for this estimation is as follows:

$$L={\displaystyle \frac{1}{2}}\sum _d\mid f_{gmm}\left(x_d\right)-f_{target}\left(x_d\right)\mid ▵x_d$$

This causes the result in Table 4 be slightly different from the results in Table 1, Table 2, and Table 3

This experiment demonstrates that the proposed method is superior to both EM and BVI. Our method, with a component size of 200, had the lowest TVD as well as the lowest training time compared to the other methods. Various factors and fundamentals can cause the EMs and BVI to underperform when the number of Gaussians increases. We consider one of the most important factor is $\sigma$. How we handle $\sigma$ is also the most notable difference between our method and others. Learning $\sigma$ throughout the training process can cause degeneration. When there are insufficiently many points for a mixture components, $\sigma$ may tend toward zero and infinite likelihood may occur or the mixture component may fall into singularity. This requires artificial regularization of the covariance and causes instability when applying the GMM in other models, such as neural networks. In the next section, we describe four applications that demonstrate how a neural network can benefit from the proposed method.

4. Neural Network Application

In this section, three neural network applications are designed to demonstrate how the GMM expansion and proposed learning algorithm can be integrated into neural networks to deliver various functionalities. The first application shows how the model performance can be enhanced when learning probabilistic outputs. The second and third applications utilize the proposed method to transform latent variables into random variables. In the second application, latent variables in an Auto-encoder [25] are converted into random variables modeled as GMM distributions. In the third application, the Gram matrices in the style transfer algorithm [9] are presented as Gram matrix distributions. In addition to the experimental results in this section, additional experiment results are given in Appendix C.

4.1. Probabilistic Output

Generally, data is not normally distributed, which makes modeling with $f\left(x\right)=y+ϵ$ inadequate. One solution is to introduce a Gaussian Mixture Model (GMM) to predict the distribution for each input, so $f\left(x\right)=GMM(\pi ,\mu ,\sigma )$. This approach addresses the issue where x and y are not unique one-to-one mappings. Tensorflow-probability2 provides a well-built package and examples to demonstrate how to implement GMM in a neural network. However, training to predict GMM parameters can be difficult. Specifically, the singularity problem may arise, as mentioned in earlier Section 3. To put it simply, $\sigma$ can approach zero during training. Probabilistic models are usually trained with likelihood functions, and when $\sigma$ tend toward zero, it can lead to infinite likelihoods, causing the cost function to explode. By applying GMM expansion, we can solve this problem and achieve better learning outcomes.

Figure 5 illustrates this concept. The classic Swiss roll data exemplify the non-unique mapping of x and y. By introducing some shifts and variance, we obtain the dataset shown on the right side of Figure 5. The upper plot shows the learning results using the standard GMM approach provided by Tensorflow-probability. In this method, the neural network is trained to output $\pi ,\mu ,\sigma$ to form a GMM for each prediction. The lower plot shows the results of our approach. By fixing $\mu$ and $\sigma$ as constants, our GMM expansion only requires a set of $\pi$ to form the GMM. The neural network in our method is trained to predict a set of $\pi$. The same likelihood function is used for training and no singularity problem should be addressed. The results are evident: the predictions from our model are visually and significantly better than those from the standard approach. In this experiment, both models are treated equally, using the same multilayer perceptron with identical hidden layer sizes and training steps. Notably, the proposed learning algorithm can also be used to preprocess data into $\pi$s. Training from the likelihood of GMM transforms into the likelihood of the discrete distribution $\pi$, which is known as entropy.

4.2. Mnist Data Set of Handwritten Digits Generation

To test whether the proposed method can be beneficial in other neural network applications, we focused on latent variables in neural networks, especially in the generative framework. In this application, the Auto-encoder [8] was tested to demonstrate the effectiveness of the proposed method. In many encoder–decoder settings, the bottleneck layer is typically responsible for controlling the generated output. The main focus of this application is to consider the latent variables for different classes as random variables, to ensure that we can learn the latent distribution with respect to different classes. Therefore, the latent distributions were treated as GMM and the proposed GMM expansion method was used to learn these distributions. Subsequently, we drew random samples from the learned GMM latent distribution to generate random images for different classes, which were digit 0-9 for the Mnist dataset. The numerical experiment was replicated from a previous study [10]. In an experiment by Kolouri et al. [10], they demonstrated how to employ the Sliced Wasserstein Distance to learn an embedding layer (latent variables) through a GMM. In this application, we learn the latent distributions using the GMM expansion and learning algorithm. The model contains an autoencoder and a straightforward convolutional classification network as the discriminator [24,25,26,27,28]. The entire process is illustrated in Figure 6.

Note that, in our training process, we slightly amended the loss function of GAN. The loss function for this application is expressed by Equation (17).

$$J^{\left(G\right)}=-\alpha {\displaystyle \frac{1}{2}}{E}_{z}log\left(1-D\left(G\left(z\right)\right)\right)+{\Vert x-G\left(z\right)\Vert }^2$$

(17)

where $J^{\left(G\right)}$ is the loss of a generator, $\alpha$ is a constant to adjust the level of discriminator loss, ${\displaystyle \frac{1}{2}}{E}_{z}log\left(1-D\left(G\left(z\right)\right)\right)$ is the discriminator loss same as [28] and x is the target image and $G\left(z\right)$ is output image generated by generator. MSE is set as the major source for controlling the direction of learning, and the GAN loss plays a less important role.

Following up the notation provided in [28], for the generator loss, we multiply $\alpha$ to the cross-entropy loss to regulate the discriminator. The addition of Mean Square Loss guides the generator to generate images close to the dataset. When the autoencoder was trained, the values of the latent variables from different classes were gathered, and the proposed method was used to learn the GMMs as latent distributions for each class. After obtaining the latent distributions, the decoder was trained independently with a mean-square error to further match each image. Subsequently, a generator based on the GMM is acquired. Images can be generated by randomly sampling the GMM from different classes. The random samples generated by our model and samples from [10] are shown in Appendix Figure A6. Contrary to the red-highlighted EM-GMM and SW-GMM results, the majority of random samples in our results did not have significant unidentifiable generations, and each class was well separated. The generations across different classes are provided in Appendix Figure A7. Two latent samples are drawn from two different classes and randomly mixed to form a new latent sample for image generation. For instance, the latent in the upper-left subplot is a combination of latent samples from 0 and 4. The output images resembled a 9. These results indicate that the latent distributions learned by the proposed method can separate generations from different classes, and a mixture of latent distributions can retain input category attributes and generate new undefined images.

4.3. Style Transfer

To further demonstrate how latent variables can be transformed into distributional random variables, we modified the style transfer algorithm [9] by converting the Gram matrix into a GMM. The Gram matrix which forms the basis of style transfer, measures the degree of correlation between extracted features as computed by a neural network. A trained VGG19 [29] feature extraction module was applied for feature extraction, as outlined in the original study. When an image is passed through VGG19, it produces output features from each layer. The Gram matrix for the features from each convolutional block was calculated using Equation (18).

$$G_{ij}^l=\sum _k{F}_{ik}^l{F}_{jk}^l,$$

(18)

A specific Gram matrix determines the specific algorithmic result. To provide diversity to the algorithm, we can change the learning rate or the weights assigned to style and content features, although performing aggressively can lead to chaotic results. Instead, a distributional Gram matrix can be used, and random sampling of the Gram matrix distributions can produce variations effortlessly. Moreover, similar to application 2, mixing latent across different classes can produce uncategorized images, The random mixture of gram matrix distribution samples from style and content can stochastically create new gram matrix to create variations.

To learn a gram matrix distribution for style and content images, images were randomly cropped and sent to VGG19 independently to obtain sets of feature values. These sets of feature values produce sets of Gram matrices, and the proposed method can be used to learn Gram matrices as GMM distributions. Finally, the same optimization procedure as in the original study [9] was performed using samples from the learned GMM gram matrix distributions. The number of style features in the Gram matrix was determined by the amount of cropped content and style images. For example, we arbitrarily cropped 25 photos from Style 1 and 75 images from Style 2. The sampled Gram matrix will probably comprise $25\%$ of the values from style 1 and $75\%$ of the values from style 2, as we sample the Gram matrix distributions. Our method is illustrated in Figure 7. Two well-known Vincent Van Gogh works Night" and "Sunflowers" are blended by this algorithm. Appendix Figure A7 shows an increasing number of style mixings in the dataset. This distributional approach changes the fine attributes of the target images but maintains the over style. The output images can be stitched together to create a short video3.

5. Conclusions, Discussion and Future Work

In this study, we present a simple yet effective method called GMM expansion for approximating arbitrary densities. Our method, along with the proposed algorithm, demonstrates competitive performance compared to state-of-the-art techniques such as the Bayesian variational Gaussian mixture method and the EM algorithm. The GMM expansion delivers better accuracy and learning speed. We have found that this method is particularly advantageous for neural networks, as illustrated through three applications. The first application demonstrates that our method is more suitable for modeling distributional outputs in neural network, which helps address non-Gaussian problems. The second and third applications demonstrate how latent variables can be effectively transformed into random variables using the proposed method.

There are several theoretical challenges and areas for further development. We have previously extended this method to high-dimensional densities, which could lead to novel embedding techniques or applications to Gaussian processes. In addition, using GMM to define latent variables as random variables relates to cardinality problems, which involve projecting latent distributions onto data distributions with a neural network. Further exploration of these theoretical problems may contribute to developing advanced generative models, such as incorporating our method into diffusion models [30,31,32].