Learning General Gaussian Mixture Model with Integral Cosine Similarity

1 Introduction

Gaussian Mixture Model (GMM) is a parametric probability density function consisting of a weighted sum of Gaussian distribution functions [Sorenson and Alspach, 1971]. It can be applied to approximate any complicated density defined on Rd with a sufficiently large number of mixture components [Kolouri et al., 2018]. GMM is a popular and powerful tool for a wide range of applications, such as pattern classification [Duda et al., 1973], statistical modeling [Zoran et al., 2012], unsupervised anomaly detection [Qu et al., 2020], etc.

Expectation Maximum (EM) algorithm is the most widely used method to estimate a GMM model [Dempster et al.,

∗Corresponding author

1977; Ng, 2000]. However, in many practical applications, the covariance matrices of some Gaussian components of the mixture model are singular (referred to as the singularity problem in the remaining text) due to insufficient data or data correlation. As the inverse of a singular matrix does not exist and the normalization term of the Gaussian distribution function is infinite, this will cause the EM algorithm to fail. A variety of methods have been proposed to overcome this difficulty. One kind of method uses data transformation as a preprocessing step and then performs modeling in the transformed space. PCA (principal component analysis) and feature selection methods are representatives of this category [Xu et al., 2014; Venegas et al., 2019]. Another kind of method, i.e., regularization, has been proposed to handle the problem by modifying the objective function or constraining the structure of covariance matrix. For instance, small values on the diagonal entries of the covariance matrix is penalized [Zong et al., 2018], or a positive quantity is added to the diagonal of the covariance matrix [Do and Ohsaki, 2021; Davari et al., 2018; Cho et al., 2019]. Another way is to restrict the covariance matrix to be diagonal [Ban et al., 2018; Ma et al., 2019; Yang et al., 2019]. Although these methods overcome the singularity problem, they may lead to the loss of information or limit the powerful representation capability of GMM [Zong et al., 2018; Zhou et al., 2020].

This paper proposes a new learning method called G2M2 (General Gaussian Mixture Model), i.e., being capable of working in general conditions even when the covariance matrices are singular. The idea is straightforward. Starting from an unnormalized Gaussian mixture function (UGMF) without a normalization term compared to GMM, we move the function towards a target but unknown data density distribution until the function and the data density distribution are at the same orientation, but may only differ by a scaling constant. Normalizing the UGMF by making its integration to 1, the unnormalized function then becomes the data density distribution we want to estimate. The key to realize this idea is the computation of the similarity between the UGMF and the unknown target data density distribution. We tackle this by introducing an integral cosine similarity (ICS) function and using Monte Carlo (MC) sampling to compute the cosine similarity without the need to explicitly estimate the data density distribution. The schematic illustration of the idea of G2M2 is shown in Fig. 1 and summarized as follows.

• Step 1 (Fig. 1(a)): starting from an UGMF ˜ p (see Eq. (3)) and assuming the unknown target probability density distribution of the data is q, we resort to the angle between ˜ p and q, i.e., θ. We introduce the ICS to compute cos2 (θ) through MC sampling without the need to explicitly estimating q (see Section 3.2).
• Step 2 (Fig 1(b)): move ˜ p towards q to maximize cos2 (θ). When cos2 (θ) = 1 or θ = 0, ˜ p and q will overlap and only differ by a scaling constant.
• Step 3 (Fig 1(c)): scale ˜ p by a constant such that its integration is 1, and we obtain p, the density distribution of the data which is in the form of a Gaussian mixture (see Section 3.3).

The contributions of this paper are summarized as follows.

• A learning method called G2M2 (general Gaussian mixture model) is proposed. It can fit an unnormalized Gaussian mixture function (UGMF) to data distribution thus overcoming the limitations associated with traditional GMM learning methods when some covariance matricies of the estimated Gaussian are singular.
• An integral cosine similarity (ICS) metric is introduced. It overcomes the limitation of Kullback-Leibler divergence or Jensen-Shannon divergence that can only evaluate the similarly between two probability density functions.
• A Siamese neural network (NN) architecture is constructed to optimize an ICS-based optimization function for learning the parameters of UGMF. With Monte Carlo sampling, the parameters of UGMF can be estimated in an elegant way.
• Experimental results show that the proposed solution outperforms conventional EM-based methods in estimating the probability density functions of complex data and achieves state-of-the-art performance in the application of unsupervised anomaly detection.

2 Preliminary and RelatedWork

Fundamentals of GMM. AGMM is defined as

where x ∈ Rd is a d-dimensional real-valued vector; λ = {ωi,µi,Σi(i = 1, · · · ,M)} is a set of parameters characterizing the GMM; ωi (i = 1, · · · ,M) is the mixture weight with the constraint ∑Mi=1 ωi = 1(ωi > 0); and g(x|µi,Σi) (i = 1, · · · ,M) is a multivariate Gaussian density component parametrized by a mean vector µi and a covariance matrix Σi, which is given as:

Note that if the matrix Σi is singular in (2), the normalization term 1 (2π) d2 |Σi| 1 2 = 10 becomes infinite and Σ-1 i does not exist. This happens due to insufficient data or data having correlations, which will result in the failure of EM.

Estimation of GMM parameters with EM. The EM algorithm [Dempster et al., 1977; Moon, 1996] is effective in estimating the GMM parameters. However, the algorithm suffers from some drawbacks, including the failure to work when the covariance matrix is singular, and being sensitive to noise and the initialization of parameters [McLachlan and Krishnan, 1997; Xu and Wunsch, 2005; Rousseeuw and Hubert, 2011]. We will show that these drawbacks can be overcome by our method.

Regularized GMM methods. The most popular way to solve the singularity problem is by adding a positive value to the diagonal of the covariance matrix [Do and Ohsaki, 2021; Davari et al., 2018; Cho et al., 2019]. However, this is equivalent to adding noise to the covariance matrix and reduces the estimation accuracy. Some other methods assume data features are independent of each other and thus setting the covariance matrices to be diagonal [Ban et al., 2018; Ma et al., 2019; Yang et al., 2019]. However, such a regularization strategy overlooks correlation in data, which is nonnegligible in many practical applications.

3 ProposedMethod

3.1 Unnormalized GaussianMixture Function

A UGMF is defined as

where η = {ξi,µi,Σ-1 i (i = 1, · · · ,M)} is a set of parameters characterizing the UGMF, and ˜ g(x|µi, Σ-1 i ) is a Gaussian function which removes the normalization term 1 (2π) d2 |Σi| 1 2 as:

Note that there will be no infinite term in ˜ p(x|η) even when the covariance matrix is singular, and we resort to direct estimation of Σ-1 i .

3.2 Integral Cosine Similarity

Let q(x) be the ground truth probability density function of data. Let Q = [q(x1), q(x2), · · · , q(xN)] and ˜ P = [˜ p(x1|η), ˜ p(x2|η), · · · , ˜ p(xN|η)] be two vectors, where xi (i ∈ {1, 2, · · · ,N}) are uniformly sampled from Ω(Rd ). As there are infinite samples (N →∞) , the squared cosine similarity of two functions ˜ p and q is defined as:

Note that unlike Kullback-Leibler divergence or JensenShannon divergence that can only evaluate the similarly between two probability density functions, Eq. (4) can evaluate the similarity between two arbitrary functions.

Define integrals I˜ p,q = ∫Ω ˜ p(x|η)q(x) dx, I˜ p2 = ∫Ω ˜ p2 (x|η) dx, and Iq2 = ∫Ω q2 (x) dx. With MC sampling, I˜ p,q can be approximated by the average of the outputs of the function ˜ p(x|η) with samples following the distribution of q(x). Consequently, given the data xri (i ∈ {1, 2, · · · ,Nr}) which are drawn from real distribution q(x), we have

I˜ p2 can be computed similarly. Let xui (i ∈ {1, 2, · · · ,Nu}) be random samples drawn from a uniform distribution defined on Ω(Rd ). For simplicity and without loss of generality, we assume each dimension in Ω is bounded, and the minimum and the maximum value among all dimensions are A and B, respectively. We adopt the MC sampling method for computing integral by inputting the samples xui into ˜ p(x|η):

In other words, I˜ p2 can be approximated by the average of the squared outputs of ˜ p(x|η) with a multiplicative constant (B -A)d . In this way, the cosine similarity can be obtained:

As q(x) and data are given, Iq2 and (B -A)d I are constant and they can be ignored. Hence we define a metric called integral cosine similarity (ICS) as:

It can be used to measure the similarity between the UGMF and the data probability density function so as to directly learn the parameter η of ˜ p(x|η) from data without the need to compute matrix inverse, which is rather different from traditional learning methods such as EM. In implementation, 1/ICS is used for optimization.

3.3 Scaling Constant

After the 1/ICS is minimized, we need a scaling constant, which is defined as C = ∫Ω ˜ p(x|η)dx, to normalize the UGMF ˜ p(x). Resorting to MC sampling again, it can be approximated by

The resultant estimated density distribution after normalization is p(x|η) = 1C ˜ p(x|η).

3.4 A Siamese Network for Learning Parameters

This paper implements all the calculations in (8) with deep learning operations so that the optimization process can be supported by existing deep learning systems. To this end, we decompose ˜ p(x|η) into some general linear and non-linear operations available in existing deep learning platforms. The decomposed form of ˜ p(x|η) can be expressed as:

where Bi = [bTi,1 , bTi,2 , · · · , bTi,d ]T is a matrix composed of d vectors, and [ti,1, ti,2, · · · , ti,d]T = Biui. The objective function (8) can be rewritten as:

Note that bTi,v xrj and bTi,v xuj are inner products, which can be implemented by convolution. ti,v is an additive parameter, and therefore it can be set as a bias term. 1 Nu ∑Nu j=1 (·) and 1 Nr ∑Nr j=1 (·) can be implemented by global average pooling (GAP). A Siamese NN structure is used and shown in Fig. 2. Note that Sβ(η) and Sγ(η) can share the same branch to obtain ˜ p(xui |η) (i ∈ {1 · · · ,Nu}). The overall learning process can be realized by minimizing the objective function Eq.(11) by using standard back propagation.

4 Experiments

In this section, the characteristics of the proposed method are demonstrated by experiments, where the Adam [Kingma and Ba, 2014] optimizer was used. The method was also applied to unsupervised anomaly detection. For characteristic evaluations, the learning rate was set to 0.1 and decayed by 0.9 every 20 epochs. For unsupervised anomaly detection, the learning rate was set to 0.01 and decayed by 0.1 every 20 epochs. The total number of training epochs was set to 40, 000 for characteristic evaluations and 80, 000 for unsupervised anomaly detection. The parameters in GMM or UGMF are randomly initialized. In this study, G2M2 was implemented 1 with Tensorflow (version 1.13.1) and run on a computer equipped with Intel Xeon E5-2640 CPU, 252 GB memory, and NVIDIA 1080Ti GPU (11 GB memory).

1 The implementation code is available at https://github.com/media-sec-lab/G2M2

4.1 Characteristic Evaluations

In this part, we first evaluate the performance of G2M2 in density estimation for 1-D data that drawn from some typical distribution functions. Then we visualize the the learning process with the proposed ICS metric for fitting 2-D Gaussian mixture. Finally, we show the robustness against noisy data and the insensitivity to parameter initialization. In the experiments, the parameters were set as: Nu = 106 , A = -30, B = 30 for 1-D data and Nu = 106 , A = -12, B = 12 for 2-D Gaussian mixture data.

UGMF and Normalized Density in Density Estimation. In our method, the parameters of UGMF ˜ p(x|η) are estimated first, and then the normalization term C is calculated to obtain p(x|η). In this way, our method does not require regularization like traditional methods. In this study, two toy examples were set up to demonstrate this characteristic. G2M2 was applied to two typical 1-D data distribution functions (with Nr = 60, 000), including Rayleigh (R) and Uniform (U). Results are shown in Fig. 3. It can be seen that ˜ p(x) is close to true distribution and only differs by a constant. When ˜ p(x) was normalized, p(x) was close to the true density with M Gaussian components.

Learning with the ICS. We propose to minimize 1/ICS (Eq. (11)) for learning the G2M2 parameters. In this experiment, we show that as the 1/ICS is minimized, the estimated density is close to the true density. We used random samples (Nr = 60, 000 and d = 2) drawn from a Gaussian mixture distribution (M = 200), in which 1% of the samples were replaced with some correlated data as noise samples, which may lead to a singular covariance matrix. When the sample number Nu approaches infinite, the ICS can well approximate squared cosine similarity. However, the sample number Nu in implementation is limited, indicating that the obtained 1/ICS has an estimation error. We may use the averaged 1/ICS as the ground truth. Figure 4 shows that the obtained 1/ICS fluctuates around the averaged value. As the averaged 1/ICS converges, the estimated density by G2M2 converges to the true density.

Robustness Analysis. The robustness of the estimation algorithm is very important. In this experiment, the robustness of the proposed method was evaluated in two aspects, including the sensitivity of ICS to noise and the sensitivity of G2M2 to parameter initialization. The similarity between the estimated density and the true density is evaluated by ALL (average log likelihood) [Kolouri et al., 2018] which is commonly employed by EM-based methods. Besides, as indicated in Section 3.2, Eq. (4) can evaluate the similarity between two arbitrary functions; therefore we can use the ICS as a kind of metric by replacing the UGMF in Eq. (8) with the density function. Without confusion, we still call it as ICS.

Robustness to Noise: To get a sense of the robustness to noise, we studied a simple scenario by using 1000 clean samples randomly drawn from N(µ,Σ), where µ = [0, 0], and Σ = [0.02 0 0 0.02]. We set M = 1 and assume Σ is known.

As a result, there is only one unknown parameter µ in Eq. (8). In this way, the 1/ICS value is a function of µ and it can be visualized. So is ALL. It can be seen from Fig. 5 that for clean data, the global optimumwas the same for both metrics. It can also be seen that when there is some noise, the contour map of ALL changes significantly, which leads to its global optimum shifts greatly from the ground-truth of µ. In contrast, the contour map of ICS was less sensitive to noise. The reason is that ALL focuses on fitting every single sample, including a noisy sample. For example, if a sample xi leads to p(xi) = 0, then log(p(xi)) = -∞, which indicates that the penalty for ALL is infinite. In this case, even the change of a few samples may cause a significant change in ALL. However, ICS pays more attention to the overall sample distribution. According to (8), only when ˜ p(xi) = 0 holds for all samples xi, the value of 1/ICS will become infinite, making it be more robust to noise.

Insensitivity to Parameter Initialization: We performed an experiment to demonstrate G2M2 is insensitive to random initialization of parameters. Our method was compared with EM-based methods, including EM-GMM-R by adding regularization values to the diagonal of the covariance matrix, and EM-GMM-D by setting the covariance matrix to be diagonal. The regularization parameter of EM-GMM-R was empirically selected from 10-6 to 10-1 and we found 10-4 achieved the best performance. In this scenario, the 2-D Gaussian mixture data polluted with noise, which was previous used and shown in the top-right corner of Fig. 4, was used again. We learned the optimized parameters from each run of random initialization using G2M2 and EM-based algorithms. We repeated 30 runs. Some estimated density functions are illustrated in Fig. 6. Moreover, we show the histograms in terms of ALL and 1/ICS to evaluate the discrepancy between the fitted model and the true data distribution. It can be seen that G2M2 achieved better performance in most runs of initialization, showing its good robustness against parameter randomization.

4.2 Unsupervised Anomaly Detection

In the application of unsupervised anomaly detection, a GMMmodel is used to learn the distribution of normal samples, while anomalies are identified as low probability samples. To evaluate the unsupervised anomaly detection performance, we conducted experiments on the datasets for outlier detection 2 , including Lympho, Cardio, Annthyroid, Shuttle, Pima, Pendigits, Satimage2, Arrthymia, Musk, Mnist, and Optdigits. Our method was compared with conventional and state-of-the-art methods, including EM-GMM-R, EM-GMMD, iForest [Liu et al., 2012], AvgKNN [Angiulli and Pizzuti, 2002], FOD [Kriegel et al., 2008], MOGAAL [Liu et al., 2020], and DAGMM [Zong et al., 2018]. The performance of these methods was evaluated by the area under curve (AUC). For both of G2M2 and EM-based GMM, the training data were normalized to [-1, 1] and the same scaling factor was used for the testing data. The parameters of G2M2 were set as follows: Nu = 106 , A = -1, B = 1, and M = 8.

2 http://odds.cs.stonybrook.edu/

Detection Performance. In the first set of experiments, the parameter setting in [Liu et al., 2020] with completely clean training data was adopted. In each run, 60% of the data taken by random sampling were used for training, and the remaining data were used for testing. It should be noted that only the data samples from the normal class were used to train the models. Table 1 shows the AUC of G2M2 and its counterparts. It can be seen that the performance of G2M2 was higher than that of the EM-based methods by 4% ∼ 5% on average and 13% ∼ 14% at most. Thus, G2M2 is more competitive than EM-based methods.

Since GMM is too complex to handle high-dimensional data, some methods such as DAGMM [Zong et al., 2018] use an autoencoder to reduce feature dimension, and optimize the GMM jointly with the autoencoder. In this study, the EM learning process in DAGMM was replaced with G2M2 . Specifically, the output of the last module of autoencoder, i.e., the estimation network in DAGMM, is the input of G2M2 , and the EM-based optimization term was replaced with (11). The adapted method was called DA-G2M2 . It can be seen in Table 1 that the performance of DA-G2M2 was significantly improved. In this case, the AUC was improved by nearly 8% on average compared to DAGMM, showing DA-G2M2 is highly competitive as compared with the-state-of-art methods.

Robustness to Noise. In the second set of experiment, the sensitivity of G2M2 to contaminated training data was investigated. In each run, 60% of the samples from the normal class were used. In addition, some samples from the anomaly class, of which the amount was equal to R% of the amount of normal training samples, were mixed into the training set as contamination noise. Fig. 7 illustrates the AUC of EMbased GMM and G2M2 . It can be seen that G2M2 achieved a stable detection accuracy when the contamination rate was increased from 1% to 5%. The results indicate traditional methods require a high quality training dataset (i.e., clean data with no contamination) to achieve good detection accuracy, while our new method G2M2 is much robust and less sensitive to the contamination of abnormal data.

5 Concluding Remarks

In this paper, a novel approach called G2M2 is proposed to overcome the singularity problem of conventional learning methods for GMM. Unlike traditional methods by transforming the data or constraining the form of the estimated covariance matrix, this paper presents a general GMM learning method that can relax the assumption of nonsingularity. Meanwhile, an optimization function based on a new metric called ICS is proposed to facilitate the parameter estimation in G2M2 . The effectiveness and robustness of G2M2 have been verified through experiments. While the new method, equipped with a Siamese neural network for learning parameters, is simple and not complicated compared to many existing deep learning architectures such as ResNet, Transformer, etc., it is computationally more demanding than EM-based solutions. In future work, we will try to overcome this problem especially for high-dimensional data.

Acknowledgments

This work was supported in part by NSFC (Grant 61872244, 62072313), Guangdong Basic and Applied Basic Research Foundation (Grant 2019B151502001), Shenzhen Program (Grant JCYJ20200109105008228), and Alibaba Innovative Research (AIR) Program.