Algorithm Learning：Clustering - Gaussian Mixture Model with Expectation-Maximization

Introduction

Gaussian mixture model（GMM）is a widely used clustering algorithm in the industry, which uses Gaussian distribution as the parametric model and utilizes maximum Expectation (EM) algorithm for training. At the same time, it’s also a method to accurately evaluate the accuracy of the posterior probability.

---

---

Induction Bias

Before describing the GMM, the induction bias needs to be explained. Whether it’s a machine or a person, the learning process can be viewed as an “induction process,” and induction requires some hypothetical preconditions. In machine learning, a learning algorithm also has a Hypothetical Preconditions, which is called induction bias.In my understanding, there is a connection and difference between induction bias and Dr. Ng’s hypothesis representation. Both are inductive premises of the machine, while induction bias is more general: it represents only a deterministic condition, not one that must be expressed as a function.

The setting of induction bias needs to be regulated. Failure to induce bias or bias is too general will result in Overfitting, while excessive restriction of induction paranoia will often result in some absurd results.The difficulty lies in finding a balance.

GMM has its very apt induction bias, it assumes that data obey the Mixture Gaussian Distribution, that is, data can be thought of as being generated from several Gaussian distributions. This bias is reasonable. It has nice properties in calculation. Moreover, by increasing the number of models, we can arbitrarily approximate any continuous PDF.

---

Theories

GMM is a simple extension of the Gaussian model. GMM utilizes a combination of multiple Gaussian distributions to characterize the data distribution. GMM is often used to solve cases where the data in the same collection contains many different distributions.The expected maximum algorithm is used to train the Gaussian mixture model.

In clustering, picking a random point from a GMM distribution can actually be divided into two steps: First, random selection is carried out among K Gaussian models, among which the probability of each Gaussian distribution being selected is actually its weight coefficient. Then only the random point selection problem in this distribution is considered, which is converted to a single Gaussian distribution problem.

The mixed Gaussian model gives you more information than the K-means model. The k-means result is that each data point is assigned to one cluster, and GMM gives the probability of assigning the data points to each cluster. This is also called a soft assignment

---

Procedure Explanation

1. Hypothesis Representation

GMM refers to a PDF model with the following form:
$$p(x;\theta )=\sum _{k=1}^K{\pi }_k\varphi (x;{\theta }_k)(1)$$where, ${\pi }_k$ is the coefficient, ${\pi }_k>=0$ and ${\sum }_{k=1}^K{\pi }_k=1$, ${\theta }_k=({\mu }_k,{\sigma }_k)$.
$\varphi (x\mid {\theta }_k)$ is Gaussian PDF:
$$\varphi (x\mid {\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}exp(-\frac{(x-{\mu }_k{)}^2}{2{\sigma }_k^2})(2)$$Each GMM consists of K Gaussian distributions, and each Gaussian is called a “Component”, these components are linearly added together to form the PDF of the model. The probability that each Component is selected is actually its coefficient ${\pi }_k$.

When GMM is used for clustering, we estimate the process of the parameter $({\pi }_k,{\mu }_k,{\sigma }_k)$, also known as “parameter estimation”, based on the existing data and assumed probability density function form.

2. The Dilemma of Maximum Likelihood Estimation

When we choose the PDF as the model, we often think of using the maximum likelihood estimation to solve the problem in the first place. But why does GMM face the “dilemma” when adopting this method?

For the training data set T = { x ( 1 ) , x ( 2 ) , . . . , x ( m ) } , x ( i ) ∈ χ ∈ R n T=\lbrace x^{(1)},x^{(2)},...,x^{(m)}\rbrace,x^{(i)}\in\chi\in R^n T={x(1),x(2),...,x(m)},x(i)∈χ∈Rn, suppose $x^{(i)}$ comes from a GMM, logarithmic likelihood function:
$$L(\theta )=\sum _{i=1}^{m}logp(x^{(i)};\mu ,\sigma ,\pi )=\sum _{i=1}^{m}log\sum _{k=1}^K{\pi }_{k}N(x^{(i)};{\mu }_k,{\sigma }_k)(3)$$
At this point, there is an accumulation sign in the logarithm sign, and it is difficult to continue the derivative calculation, so maximum likelihood estimation should be in trouble.

3. The Bayesian Understanding of GMM

As we mentioned, ${\pi }_k$ can be viewed as the probability that the $k$th class is selected. We introduce a $K$ dimensional random variable $z=(z_1,z_2,...,z_K)$, $z_k(1\le k\le K)$, $z_k=1$ indicates that the $k$th gaussian distribution is selected, and $z_k=0$ means that it is not. It has the following properties:
z k = { 0 , 1 } z_k=\lbrace 0,1\rbrace zk​={0,1}$$p(z_k=1)={\pi }_k$$$$\sum _K{z}_k=1$$We assume that $z_k$ are independently identically distributed, we can write the joint probability distribution form of $z$:
$$p(z)=p(z_1)p(z_2)...p(z_K)=\prod _{k=1}^{K}p(z_k)=\prod _{k=1}^K{\pi }_k^{z_k}(4)$$Where $z$ can only be 0 or 1, and only $z_k$ is 1, while all the other $z_j(j\ne k)$ are 0.

Each class of data selected is Normally Distributed, which can be further represented by conditional probability:
$$p(x\mid z_k=1)=N(x;{\mu }_k,{\sigma }_k)(5)$$$$p(x\mid z)=\prod _{k=1}^{K}N(x;{\mu }_k,{\sigma }_k{)}^{z_k}(6)$$According to (4) and (6), combined with conditional probability, we can get the form of $p(x)$:
$$p(x)=\sum _{z}p(x,z)=\sum _{z}p(z)p(x\mid z)$$$$=\sum _z\prod _{k=1}^K({\pi }_k^{z_k}N(x;{\mu }_k,{\sigma }_k{)}^{z_k})(7)$$$$=\sum _{k=1}^K{\pi }_{k}N(x;{\mu }_k,{\sigma }_k)(8)$$It can be found that (7) has the same form as the original mixed Gaussian model (8), and the latent variable $z$ is introduced into (7).

latent variable: We know that the data can be divided into two classes, but if we randomly select a point, we don't know whether the data point belongs to the first class or the second class, and we can't observe its attribution, so we introduce latent variable to describe this phenomenon.

Based on Bayes’ idea, we can find the posterior probability $p(z\mid x;\theta )$:

$$p(z\mid x;\theta )=p(z_k=1\mid x)$$$$=\frac{p(z_k=1)p(x\mid z_k=1)}{p(x)}$$$$=\frac{p(z_k=1)p(x\mid z_k=1)}{{\sum }_{i=1}^{K}p(z_i=1)p(x\mid z_i=1)}$$$$=\frac{{\pi }_{k}N(x_k\mid {\mu }_k,{\sigma }_k)}{{\sum }_{i=1}^K{\pi }_{i}N(x;{\mu }_i,{\sigma }_i)}(9)$$We define the above formula (9) $\gamma (z_k)$ to represent the posterior probability of the $k$th Component and lay a foundation for the following Expectation-Maximum algorithm.

4. Expectation-Maximum（EM） Algorithm

EM algorithm is an iterative algorithm for maximum likelihood estimation or maximum posterior probability estimation of probabilistic model parameters with latent variables. Each iteration of EM algorithm is accomplished in two steps: Step E, expectation; Step M, maximization.

EM algorithm optimizes the parameter step by step by iterating to maximize $L(\theta )$. A “clever” method is to find the lower bound of $L(\theta )$. As long as the lower bound of $L(\theta )$ is maximized during each iteration, $L(\theta )$ can be guaranteed to increase all the time.

Jensen Inequality

In order to solve the above problems, the Jensen Inequality is introduced:
$$f(\sum _{i=1}^N{\lambda }_i{x}_i)\le \sum _{i=1}^N{\lambda }_{i}f(x_i)(10)$$$$s.t.{\lambda }_i\ge 0,\sum _{i=1}^N{\lambda }_i=1$$Where $f(x)$ is convex function. If and only if $Xi\equiv C$(constant $C$) the equal sign is true.

And the probability form of the inequality is shown as follows:$$f(E(X))\le E(f(X))(11)$$

Step E Algorithm Derivation

We utilize Jensen’s inequality to solve the lower bound of $L(\theta )$ in the $t+1$ iteration. The scaling process is as follows:
$$L(\theta )\ge \sum _{i=1}^m\sum _{z^{(i)}}T(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{T(z^{(i)})}(12)$$
Where $T(z^{(i)})$ is a constructed ${\lambda }_i$ in (10), it must be a known when the estimator ${\theta }^{(i)}$'s quantity is given.

According to the Bayesian understanding above, we find that the posterior probability of $z$ exactly meets the requirements of the construction $T(z^{(i)})$ quantity, and because:$$\sum _{z^{(i)}}p(z^{(i)}\mid x;{\theta }^{(t)})=1(13)$$$$\frac{p(x^{(i)},z^{(i)};{\theta }^{(t)})}{p(z^{(i)}\mid x;{\theta }^{(t)})}=\frac{p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})p(x^{(i)};{\theta }^{(t)})}{p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})}$$$$=p(x^{(i)};{\theta }^{(t)})$$$$\equiv C(14)$$
In fact, the EM algorithm also utilize $p(z)$ to construct the lower bound function. We can construct as follows:$$T(z^{(i)})=p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})$$$$=\frac{p(z_k^{(i)}=1)p(x^{(i)}\mid z_k^{(i)}=1)}{{\sum }_{k=1}^{K}p(z_k^{(i)}=1)p(x^{(i)}\mid z_k^{(i)}=1)}$$$$=\frac{{\pi }_{k}N(x^{(i)}\mid {\mu }_k,{\sigma }_k)}{{\sum }_{k=1}^K{\pi }_{k}N(x^{(i)};{\mu }_k,{\sigma }_k)}(15)$$
$$w_k^i=p(z_k^{(i)}=1\mid x^{(i)};{\theta }^{(t)})(16)$$$w_j^i$ represents the probability that the $i$th sample comes from the $j$th Gaussian distribution, that is, $w_j^i$ is a matrix of $n$_samples $\ast$ $k$_components.

And that’s where we get the key $Q$ function（Doctor Li Hang expresses $Q$ function in Statistical Learning Method: it is the expectation of logarithmic likelihood function $p(x^{(i)},z^{(i)};\theta )$ of complete data on conditional probability distribution $p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})$ of latent data $z$ under given observation data $x$ and current parameter ${\theta }^{(t)}$）:
$$Q(\theta ,{\theta }^{(t)})=\sum _{i=1}^m\sum _{z^{(i)}}T(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta )}{T(z^{(i)})}$$$$=\sum _{i=1}^m\sum _{z^{(i)}}p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})log\frac{p(x^{(i)},z^{(i)};\theta )}{p(z^{(i)}\mid x^{(i)};{\theta }^{(t)})}(17)$$

Step M Algorithm Derivation

Once $Q$ function is obtained, the function is maximized to obtain the model parameters for the next iteration:
$${\theta }^{(t+1)}=argmaxQ(\theta ,{\theta }^{(t)})$$First, we substitute (4), (5) and (16) into the $Q$ function to get:$$Q(\theta ,{\theta }^{(t)})=\sum _{i=1}^m\sum _{k=1}^{K}T(z_k^{(i)}=1)log\frac{p(x^{(i)}\mid z_k^{(i)}=1;\mu ,\sigma )p(z_k^{(i)}=1;\pi )}{T(z_k^{(i)}=1)}$$$$=\sum _{i=1}^m\sum _{k=1}^K{w}_k^{i}log\frac{{\pi }_k\frac{1}{(2\pi {)}^{\frac{n}{2}}{\sigma }_k^{\frac{1}{2}}}exp(-\frac{1}{2}(x^{(i)}-{\mu }_k{)}^T{\sigma }_k^{-1}(x^{(i)}-{\mu }_k))}{w_k^i}(18)$$
Take the partial derivative of the above expression with respect to each parameter, set the partial derivative equal to 0, obtain the parameter value of the maximum value of the function and assign it to the parameters participating in the next iteration:
$${\mu }_k:=\frac{{\sum }_{i=1}^m{w}_k^i{x}^{(i)}}{{\sum }_{i=1}^m{w}_k^i}(19)$$$${\sigma }_k:=\frac{{\sum }_{i=1}^m{w}_k^i(x^{(i)}-{\mu }_k{)}^T(x^{(i)}-{\mu }_k)}{{\sum }_{i=1}^m{w}_k^i}(20)$$$${\pi }_k:=\frac{{\sum }_{i=1}^m{w}_k^i}{m}(21)$$Where：$$w_k^i=p(z_k^{(i)}=1\mid x^{(i)};{\theta }^{(t)})$$

EM Algorithm Flow

Input: Observational dataset T = { x ( 1 ) , x ( 2 ) , . . . , x ( m ) } , x ( i ) ∈ χ ∈ R n T=\lbrace x^{(1)},x^{(2)},...,x^{(m)}\rbrace,x^{(i)}\in\chi\in R^n T={x(1),x(2),...,x(m)},x(i)∈χ∈Rn;    \qquad\,\, Gaussian Mixture Model $p(x;\theta )$

Output: GMM parameters $(\pi ,\mu ,\sigma )\in R^K$

1. Define the number of components $K$, initialize the parameters$(\pi ,\mu ,\sigma )$, and start the iteration.
2. Step E: According to the model parameters of the current iteration $t$,calculate the responsiveness $w_k^i$ of the sub-model $k$ to the observed data $x^{(i)}$:$$w_k^i=\frac{{\pi }_{k}N(x^{(i)}\mid {\mu }_k,{\sigma }_k)}{{\sum }_{k=1}^K{\pi }_{k}N(x^{(i)};{\mu }_k,{\sigma }_k)},i=1,2,...,m;k=1,2,...,K$$
3. Step M: Calculate the model parameters of the new iteration $t+1$：$${\mu }_k^{new}:=\frac{{\sum }_{i=1}^m{w}_k^i{x}^{(i)}}{{\sum }_{i=1}^m{w}_k^i}$$$${\sigma }_k^{new}:=\frac{{\sum }_{i=1}^m{w}_k^i(x^{(i)}-{\mu }_k{)}^T(x^{(i)}-{\mu }_k)}{{\sum }_{i=1}^m{w}_k^i}$$$${\pi }_k^{new}:=\frac{{\sum }_{i=1}^m{w}_k^i}{m}$$
4. Calculate the maximum likelihood function $L(\theta )$, test whether the function converges, if not, return step 2.

---

Convergence Verification of EM Algorithm

In fact, the EM algorithm maximizes the logarithmic likelihood function by iterative stepwise approximation. Since the likelihood function exactly has an upper bound, we expect parameter iteration can increase $L(\theta )$:
$$L({\theta }^{(t+1)})-L({\theta }^{(t)})=L({\theta }^{(t+1)})-Q({\theta }^{(t)},{\theta }^{(t)})$$$$\ge Q({\theta }^{(t+1)},{\theta }^{(t)})-Q({\theta }^{(t)},{\theta }^{(t)})$$As the value of the new iteration parameter makes $Q$ function reach the maximum value, there are:$$Q({\theta }^{(t+1)},{\theta }^{(t)})-Q({\theta }^{(t)},{\theta }^{(t)})\ge 0$$So the likelihood function can converge during iteration, and EM algorithm works.

---

Potential Applications

Since normal distribution is common in daily life, GMM has high application value. Combined with big data technology, GMM can locate information about individuals’ behavior within population. For example, using facial data to identify a person’s race.

As a classification model, GMM can try to be a part of GAN. GMM improves the performance of the final generator by providing more specific categorization information.