This note represents my attempt at explaining the EM algorithm (Hartley, 1958;
Dempster et al., 1977; McLachlan and Krishnan, 1997). This is just a slight
variation on Tom Minka’s tutorial (Minka, 1998), perhaps a little easier (or perhaps
not). It includes a graphical example to provide some intuition.
1 Intuitive Explanation of EM
EM is an iterative optimization method to estimate some unknown parameters, given measurement data U. However, we are not given some “hidden” nuisance variables J, which need to be integrated out. In particular, we want to maximize the posterior probability of the parameters given the data U, marginalizing over J:
= argmax
X
J2Jn
P(; JjU) (1)
The intuition behind EM is an old one: alternate between estimating the unknowns and the hidden variables J. This idea has been around for a long time. However, instead of finding the best J 2 J given an estimateat each iteration, EM computes a distribution over the space J . One of the earliest papers on EM is (Hartley, 1958), but the seminal reference that formalized EM and provided a proof of convergence is the “DLR” paper by Dempster, Laird, and Rubin (Dempster et al., 1977). A recent book devoted entirely to EM and applications is (McLachlan and Krishnan, 1997), whereas (Tanner, 1996) is another popular and very useful reference. One of the most insightful explanations of EM, that provides a deeper understanding of its operation than the intuition of alternating between variables, is in terms of lowerbound maximization (Neal and Hinton, 1998; Minka, 1998). In this derivation, the E-step can be interpreted as constructing a local lower-bound to the posterior distribution, whereas the M-step optimizes the bound, thereby improving the estimate for the unknowns. This is demonstrated below for a simple example.
1
−5 −4 −3 −2 −1 0 1 2 3 4 5
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Figure 1: EM example: Mixture components and data. The data consists of three samples drawn from each mixture component, shown above as circles and triangles.
The means of the mixture components are 2 and 2, respectively.
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
q2
Figure 2: The true likelihood function of the two component means 1 and 2, given the data in Figure 1.
2
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
i=1, Q=−3.279564
q2
1 2 3 4 5 6
0.5
1
1.5
2
2.5
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
i=2, Q=−2.788156
q2
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
i=3, Q=−1.501116
q2
1 2 3 4 5 6
0.5
1
1.5
2
2.5
1 2 3 4 5 6
0.5
1
1.5
2
2.5
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
i=4, Q=−0.817491
q2
−3
−2
−1
0
1
2
3
−3
−2
−1
0
1
2
3
0
0.1
0.2
0.3
0.4
0.5
q1
i=5, Q=−0.762661
q2
1 2 3 4 5 6
0.5
1
1.5
2
2.5
1 2 3 4 5 6
0.5
1
1.5
2
2.5
Figure 3: Lower Bounds
3
Consider the mixture estimation problem shown in Figure 1, where the goal is to estimate the two component means 1 and 2 given 6 samples drawn from the mixture, but without knowing from which mixture each sample was drawn. The state space is two-dimensional, and the true likelihood function is shown in Figure 2. Note that there are two modes, located respectively at (2; 2) and (2;2). This makes perfect sense, as we can switch the mixture components without affecting the quality of the solution.
Note also that the true likelihood is computed by integrating over all possible data associations, and hence we can find a maximum likelihood solution without solving a correspondence problem. However, even for only 6 samples, this requires summing over the space of 64 possible data-associations.
EM proceeds as follows in this example. In the E-step, a “soft” assignment is computed that assigns a posterior probability to each possible association of each individual sample.
In the current example, there are 2 mixtures and 6 samples, so the computed probabilities can be represented in a 2 6 table. Given these probabilities, EM computes a tight lower bound to the true likelihood function of Figure 2. The bound is constructed such that it touches the likelihood function at the current estimate, and it is only close to the true likelihood in the neighborhood of this estimate. The bound and its corresponding probability table is computed in each iteration, as shown in Figure 3. In this case, EM was run for 5 iterations. In the M-step, the lower bound is maximized (shown by a black asterisk in the figure), and the corresponding new estimate (1; 2) is guaranteed to lie closer to the location of the nearest local maximum of the likelihood. Each next bound is an increasingly better approximation to the mode of the likelihood, until at convergence the bound touches the likelihood at the local maximum, and progress can no longer be made. This is shown in the last panel of Figure 3.
2 EM as Lower Bound Maximization
EM can be derived in many different ways, one of the most insightful being in terms of lower bound maximization (Neal and Hinton, 1998; Minka, 1998), as illustrated with the example from Section 1. In this section, we derive the EM algorithm on that basis, closely following (Minka, 1998). The goal is to maximize the posterior probability (1) of the parameters given the data U, in the presence of hidden data J. Equivalently, we can maximize the logarithm of the joint distribution (which is proportional to the posterior):
= argmax
log P(U;) = argmax
log
X
J2Jn
P(U; J;) (2)
The idea behind EM is to start with a guess t for the parameters , compute an easily computed lower bound B(;t) to the function log P(jU), and maximize that bound instead. If iterated, this procedure will converge to a local maximizer of the objective function, provided the bound improves at each iteration. To motivate this, note that the key problem with maximizing (2) is that it involves the logarithm of a (big) sum, which is difficult to deal with. Fortunately, we can construct
4
a tractable lower bound B(;t) that instead contains a sum of logarithms. To derive the bound, first trivially rewrite log P(U;) as
log P(U;) = log
X
J2Jn
P(U; J;) = log
X
J2Jn
f t(J)
P(U; J;)
f t(J)
where f t(J) is an arbitrary probability distribution over the space J n of hidden variables J. By Jensen’s inequality, we have
B(;t)
=
X
J2Jn
f t(J) log
P(U; J;)
f t(J)
<= log
X
J2Jn
f t(J)
P(U; J;)
f t(J)
Note that we have transformed a log of sums into a sum of logs, which was the prime motivation.
2.1 Finding an Optimal Bound
EM goes one step further and tries to find the best bound, defined as the boundB(;t) that touches the objective function log P(U;) at the current guess t. Intuitively, finding the best bound at each iteration will guarantee that we obtain an improved estimate t+1 when we locally maximize the bound with respect to . Since we know B(;t) to be a lower bound, the optimal bound at t can be found by maximizing
B(t;t) =
X
J2Jn
f t(J) log
P(U; J;t)
f t(J)
(3)
with respect to the distribution f t(J). Introducing a Lagrange multiplier to enforce the constraint
P
J2Jn f t(J) = 1, the objective becomes
G(f t) =
"
1
X
J2Jn
f t(J)
#
+
X
J2Jn
f t(J) log P(U; J;t)
X
J2Jn
f t(J) log f t(J)
Taking the derivative
@G
@f t(J)
= + log P(U; J;t) log f t(J) 1
and solving for f t(J) we obtain
f t(J) =
P(U; J;t P )
J2Jn P(U; J;t)
= P(JjU;t)
By examining the value of the resulting optimal bound at t we see that it indeed touches the objective function:
B(t;t) =
X
J2Jn
P(JjU;t) log
P(U; J;t)
P(JjU;t)
= log P(U;t)
5
2.2 Maximizing The Bound
In order to maximize B(;t) with respect to , note that we can write it as
B(;t)
=
hlog P(U; J;)i + H
= hlog P(U; Jj)i + log P() + H
= Qt() + log P() + H
where h:i denotes the expectation with respect to f t(J)
=
P(JjU;t), and
Qt() is the expected complete log-likelihood, defined as:
Qt()
=
hlog P(U; Jj)i
P() is the prior on the parameters
H
=
hlog f t(J)i is the entropy of the distribution f t(J)
Since H does not depend on , we can maximize the bound with respect to using the first two terms only:
t+1 = argmax
B(;t) = argmax
Qt() + log P()
(4)
2.3 The EM Algorithm
At each iteration, the EM algorithm first finds an optimal lower boundB(;t) at the current guess t (equation 3), and then maximizes this bound to obtain an improved estimatet+1 (equation 4). Because the bound is expressed as an expectation, the first step is called the “expectation-step” or E-step, whereas the second step is called the “maximization-step” or M-step. The EM algorithm can thus be conveniently summarized as:
E-step: calculate f t(J)
=
P(JjU;t)
M-step: t+1 = argmax [Qt() + log P()]
It is important to remember that Qt() is calculated in the E-step by evaluating f t(J) using the current guess t (hence the superscript t), whereas in the M-step we are optimizing Qt() with respect to the free variableto obtain the new estimatet+1. It can be proved that the EM algorithm converges to a local maximum of log P(U;), and thus equivalently maximizes the log-posterior log P(jU) (Dempster et al., 1977; McLachlan and Krishnan, 1997).
6
A Relation to the Expected Log-Posterior
Note that we have chosen to define Qt() as the expected log-likelihood as in (Dempster et al., 1977; McLachlan and Krishnan, 1997), i.e.,
Qt()
=
hlog P(U; Jj)i
An alternative route is to compute the expected log-posterior (Tanner, 1996):
Here the second term does not depend on J and can be taken out of the expectation, and the last term does not depend on . Hence, maximizing (5) with respect to is equivalent to (4):
argmax
hlog P(jU; J)i = argmax
[hlog P(U; Jj)i + log P()]
= argmax
Qt() + log P()
This is of course identical to (4).
References
[1] Dempster, A., Laird, N., and Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39(1):1–38.
[2] Hartley, H. (1958). Maximum likelihood estimation from incomplete data. Biometrics, 14:174–194.
[3] McLachlan, G. and Krishnan, T. (1997). The EM algorithm and extensions. Wiley series in probability and statistics. John Wiley & Sons.
[4] Minka, T. (1998). Expectation-Maximization as lower bound maximization. Tutorial published on the web at http://www-white.media.mit.edu/ tpminka/papers/em.html.
[5] Neal, R. and Hinton, G. (1998). A view of the EM algorithm that justifies incremental, sparse, and other variants. In Jordan, M., editor, Learning in Graphical Models. Kluwer Academic Press.
[6] Tanner, M. (1996). Tools for Statistical Inference. Springer Verlag, New York. Third Edition.
7
Intuitive Explanation of EM
We want to maximize the posterior probability of the parameters Θ = {θ1, …, θn} given the data U, marginalizing over J:
(1)
Instead of finding the best J ∈ J given an estimate Θ at each iteration, EM computes a distribution over the n-dimensional space Jn.
EM as Lower Bound Maximization
EM can be derived in many different ways, one of the most insightful being in terms of lower bound maximization (Neal and Hinton, 1998; Minka, 1998).
The goal is to maximize the posterior probability (1) of the parameters Θ given the data U, in the presence of hidden data J. Equivalently, we can maximize the logarithm of the joint distribution (which is proportional to the posterior):
(2)
Note that the key problem with maximizing (2) is that it involves the logarithm of a (big) sum, which is difficult to deal with.
To derive the bound, first trivially rewrite log P(U,Θ) as
where ft(J) is an arbitrary probability distribution over the space Jn of hidden variable J.
By Jensen's inequality, we have
Finding an Optimal Bound
EM tries to find the best bound, defined as the bound B(Θ;Θt) that touches the objective function log P(U,Θ).
Introducing a Lagrange multiplier λ to enforce the constraint
the objective becomes
Taking the derivative
We see that the resulting optimal bound at Θt indeed touches the objective function:
Maximizing the Bound
(3)
where $<.>$ denotes the expectation with respect to $f^t(\textbf{J}) = P(\textbf{J}|\textbf{U},\Theta^t)$.
$P(\Theta)$ is the prior on the parameters $\Theta$
$H := -<\log f^t(\textbf{J})>$ is the entropy of the distribution $f^t(\textbf{J})$
Since $H$ does not depend on $\Theta$, we consider the first two terms only:
For your reference, on the third line of (3), the last term in $<>$ does not depend on $\Theta$.
So, $<\log P(\textbf{U},\textbf{J},\Theta)>$ can be rewritten as $<\log P(\textbf{U},\textbf{J}|\Theta)> + \log P(\Theta)$.
The EM algorithm
The bound is expressed as an expectation.
Expectation step: first find an optimal lower bound B(Θ;Θt) at the current guess Θt
Maximization step: and then maximize this bound to obtain an improved estimate Θt+1.
The EM algorithm can thus be conveniently summarized as:
The Expectation Maximization Algorithm - Frank Dellaert (2002).p