 Introduction
 Traditional Derivation of EM
 An Entropy Interpretation
 DoublyStochastic EM
 Concluding Remarks
3.2 Traditional Derivation of EM
Each EM iteration is composed of two steps—Estimation (E) and Maximization (M). The Mstep maximizes a likelihood function that is further refined in each iteration by the Estep. This section derives the traditional EM and establishes its convergence property.
3.2.1 General Analysis
The following notations are adopted.

X = {x_{t} ∊ ℜ^{D}; t = 1, . . , T} is the observation sequence, where T is the number of observations and D is the dimensionality of x_{t}.

C = {C^{(1)}, . . , C^{(}^{J}^{)}} is the set of cluster mixture labels, where J is the number of mixture components.

Z = {z_{t} ∊ C; t = 1, . . , T} is the set of missing data (specifying the hiddenstate information).

θ = {θ^{(j)}; j = 1, . . , J} is the set of unknown parameters that define the density function for approximating the true probability density of X.

θ^{(j)}= {π^{(j}^{)}, φ^{(j)}}, where π^{(}^{j}^{)} denotes the prior probability of the jth component density and φ^{(j)} defines the jth component density.
Note that the combination of observations X and the "hiddenstates" Z constitute the completedata. The likelihood of the completedata is instrumental in accordance with the EM formulation.
To facilitate the derivation, define
as the loglikelihood of the incompletedata given the current estimate θ_{n}, where n represents the iteration index; also, define p(Z, Xθ_{n}) as the completed data likelihood. According to probability theory, ^{ [1] } p(Xθ_{n}) can be expressed as
Using Eq. 3.2.1 and Eq. 3.2.2, one can write the incompletedata loglikelihood as follows:
where E_{Z}{} denotes expectation with respect to Z. Thus, denote
and
where R(θθ_{n}) is an entropy term representing the difference between the incompletedata likelihood and the expectation of the completeddata likelihood. Interpretation of R(θθ_{n}) and its role in the EM algorithm is discussed further in Section 3.3.
3.2.2 Convergence Property of EM
The following demonstrates why the EM algorithm has a general convergence property. The basic idea is via Jensen's inequality. More precisely, it can be shown that if the Qfunction in Eq. 3.2.4 is improved in each iteration (in the Mstep), then so will be the likelihood function L.
The proof of convergence begins with the observation of the following relationship:
Using Eq. 3.2.6 and Jensen's inequality, this is obtained:
In the Mstep of the nth iteration, 0* is selected according to
This means one can always choose a θ* at iteration n such that
Note that this equation constitutes a sufficient condition to ensure the convergence property of the EM algorithm because, according to Eqs. 3.2.3, 3.2.7, and 3.2.9
Instead of directly maximizing L(Xθ), the EM algorithm divides the optimization problem into two subproblems: Expectation and Maximization.
In each EM iteration, the Estep computes Q(θθ_{n}) using a set of presumed model parameters θ_{n}. The Mstep determines the value of 0 (say θ^{*}) that maximizes Q(θθ_{n}); that is,
This results in (see Problem 8)
Dividing the optimization into two interdependent steps is most useful if optimizing Q(θθ_{n}) is simpler than that of L(Xθ_{n}). Figure 3.4 illustrates how the Eand Msteps interplay to obtain a maximumlikelihood solution. The next section explains how to compute Q(θθ_{n}) in the Estep and how to maximize Q(θθ_{n}) in the Mstep.
Figure 3.4 The flow of the EM algorithm.
Generalized EM
In case θ* in Eq. 3.2.8 is difficult to attain, the EM approach is still applicable if one can improve Q(θθ_{n}) in each Mstep (e.g., by gradient ascent). The algorithm is known as generalized EM. Although convergence of generalized EM is slower than that of the standard EM, it offers a more general and flexible framework for dividing the optimization process into the EM steps.
3.2.3 CompleteData Likelihood
EM begins with an optimization of a likelihood function, which may be considerably simplified if a set of "missing" or "hidden" data is assumed to be known. The following demonstrates that computing the expectation of the completedata likelihood in the Estep can be accomplished by finding the expectation of the missing or hidden data.
If X = {x_{t}; t = 1,..., T} contains T statistically independent vectors and Z = {z_{t} ∊ C; t = 1,..., T}, where z_{t} = C^{(j)} means that the jth mixture generates x_{t}, then one can write p(Z, Xθ) as
Now, a set of indicator variables is introduced to indicate the status of the hiddenstates: ^{ [2] }
where
Since for each t only one of the terms in is equal to one and all of the others are equal to 0, one can express p(Z, Xθ) as follows:
Hence, the completeddata likelihood is given by
where π^{(}^{j}^{)} is the mixing coefficient of the jth mixture. Eq. 3.2.12 uses the fact that and . Moreover, because there is only one nonzero term inside the summation , one can extract from the log function without affecting the result.
EStep
Taking the expectations of Eq. 3.2.12 and using the defintions in Eq. 3.2.4, one obtains
Then, define
and denote as the jth mixture coefficient at iteration n. Using the Bayes theorem, one can express as
The Estep determines the best guess of the membership function . Once the probability are computed for each t and j, Q(θθ_{n}) can be considered as a function of θ. In the Mstep of each iteration, this function is maximized to obtain the best value of θ (denoted as θ*). In most cases, the Mstep is substantially simplified if are known. Therefore, the Estep can be viewed as a preparation step for the Mstep.
3.2.4 EM for GMMs
To better illustrate the EM steps, a simple example applying EM to Gaussian mixture models (GMMs) is presented next. The most common forms for the mixture density are the radial basis functions (RBFs) or the more general elliptical basis functions (EBFs). In the latter case, the component density is a Gaussian distribution, with the model parameter of the jth cluster φ^{(}^{j}^{)} = {μ^{(}^{j}^{)}, Σ^{(j)}}I consisting of a mean vector and a fullrank covariance matrix.
Assume a Gaussian mixture model:
where π^{(}^{j}^{)}, μ^{(}^{j}^{)} and Σ^{(}^{j}^{)} denote, respectively, the mixture coefficient, mean vector, covariance matrix of the jth component density. The GMM's output is given by
where
is the jth Gaussian density of the GMM. A closer look at Eqs. 3.2.15 and 3.2.16 reveals that the GMM parameters θ can be divided into two groups: one containing π^{(j)}s and another containing μ^{(}^{j}^{)}s and ∑^{(}^{j}^{)}s. The former indicates the importance of individual mixture densities via the prior probabilities π^{(}^{j}^{)}s, whereas the latter is commonly regarded as the kernel parameter defining the form of the mixture density. Unlike other optimization techniques (e.g., gradient descent) in which unknown parameters can be arranged in any order, the EM approach effectively makes use of the structural relationship among the unknown parameters to simplify the optimization process.
After the initialization of θ_{0}, the EM iteration is as follows:

Estep. In the nth iteration, compute (x_{t}) for each j and t using Eqs. 3.2.14 and 3.2.16. This is followed by the Mstep described next.

Mstep. Maximize Q(θθ_{n}) with respect to θ to find θ*. Replace θ_{n} by θ*. Then, increment n by 1 and repeat the Estep until convergence.
To determine μ^{(}^{k}^{)}*, set , which gives
To determine Σ^{(k)}*, set , which gives
To determine π^{(k)}*, maximize Q(θθ_{n}) with respect to π^{(}^{k}^{)} subject to the constraint , which gives
The detailed derivations of Eq. 3.2.17 to Eq. 3.2.19 are as follows:
To determine Σ^{(}^{k}^{)}*, k = 1,..., J, let Λ(k)) = (Σ^{(}^{k}^{)})^{−1} and set , that is,
Note that Eq. 3.2.21 makes use of the identity , where A is a symmetric matrix. Note also that one can replace μ^{(}^{k)} by μ^{(}^{k}^{)}* in Eq. 3.2.20 to obtain Eq. 3.2.18.
To determine π^{(}^{r}^{)}, r = 1,..., J, maximize Q(θθ_{n}) with respect to π^{(}^{r}^{)} subject to the constraint . More specifically, maximize the function f(λ, π^{(j)}) = where A is the Lagrange multiplier. Setting results in
Summing both size of Eq. 3.2.24 from r = 1 to J, one has
Substituting Eq. 3.2.26 into Eq. 3.2.24 results in
Complexity of EM
Let T denote the number of patterns, J the number of mixtures, and D the feature dimension, then the following is a rough estimation of the computation complexity of using EM to train a GMM:
Numerical Example 1
This example uses the data in Figure 3.3(a) as the observed data. Assume that when EM begins, n = 0 and
Therefore, one has
and
Substituting X = {1, 2, 3, 4, 6, 7, 8} into Eqs. 3.2.28 and 3.2.29, Table 3.1 is obtained. Substituting in Table 3.1 into Eqs. 3.2.17 through 3.2.19 results in
Table 3.1. Values of in Example 1
Pattern Index (t) 
Pattern (x_{t}) 


1 
1 
1 
0 
2 
2 
1 
0 
3 
3 
1 
0 
4 
4 
1 
0 
5 
6 
0 
1 
6 
7 
0 
1 
7 
8 
0 
1 
Then, continue the algorithm by computing Q(θθ_{n}) that is, —which are then substituted into Eqs. 3.2.17 through 3.2.19 to obtain θ_{2.} Figure 3.5 depicts the movement of the component density functions specified by μ^{(j}^{)} and σ^{(j}^{)} during the EM iterations, and Table 3.2 lists the numerical values of Q(θθ_{n}) and θ_{n} for the first five iterations. It is obvious that the algorithm converges quickly in this example.
Figure 3.5 Movement of the component density function specified by μ(j) and (σ(j))2 for the first two EM iterations.
Table 3.2. Values of Q(θθ_{n}), μ^{(j)} and (σ^{(}^{j}^{)})^{2} in the course of EM iterations. Data shown in Figure 3.3(a) were used as the observed data.
Iteration (n) 
Q(θθ_{n}) 


0 
∞ 
0 
1 
9 
1 
1 
43.71 
2.50 
1.25 
6.99 
0.70 
2 
25.11 
2.51 
1.29 
7.00 
0.68 
3 
25.11 
2.51 
1.30 
7.00 
0.67 
4 
25.10 
2.52 
1.30 
7.00 
0.67 
5 
25.10 
2.52 
1.30 
7.00 
0.67 