EM Algorithm as I Know It
Dedicated to Prof. Dr Michail Ivanovich Schlesinger
Prof. Dr Schlesinger invented the famous unsupervised learning algorithm —the EM algorithm. More precisely, the first complete description of a discrete case of the EM algorithm known to me was published as the article "The interaction of learning and self-organization in pattern recognition" in 1968 — almost ten years before "Maximum Likelihood from Incomplete Data via the EM Algorithm" by Arthur Dempster, Nan Laird, and Donald Rubin. Michail Schlesinger did not name the algorithm “Expectation-Maximisation” — he used the term “self-organization” instead.
During his lectures, Prof. Dr Schlesinger was saying, "You should know this algorithm as if you have invented it!" While I did not reach that level of understanding of the applied part of the EM algorithm, I internalised the proof in slightly different terms than he presented in his original paper or in his lectures, where I had the honour of being a student for years (and eventually becoming a lecturer myself).
The entire algorithm inference and monotonicity proof is not too long and can be written as a single chain of equalities with one inequality within (and a couple of redefinitions for shortness):
You may start by jumping straight to The EM Algorithm because this story starts from the “reinvention” of this algorithm and proving its main quality.
Table of Contents
- Expanding Log-Likelihood
- Improving The Log-Likelihood
- Prior Probabilities
- Conditional Probabilities
- Posterior Probabilities
- Monotonicity of The EM Algorithm - Correctness of The EM Algorithm
- The EM Algorithm
- Input
- Initialisation
- Expectation (The E-Step)
- Maximisation (The M-Step)
- Result - Conclusions and Further Directions
Expanding Log-Likelihood
Our aim is, given a finite non-empty random sample of observations, to cluster it without providing any information about the desired assignments for any items of the sample.
That is why the algorithm inference starts with a random sample — a tuple of mutually independent random variables. As we know, the probability of a tuple of independent random variables is a product of their probabilities. Thus, the logarithm is a sum of logarithms of the probabilities, also known as a log-likelihood. The superscript index t means the step number because we are exploring an iterative algorithm:
What is next? We need an equivalent representation of the observation probability that preserves the equalities but moves us forward. For this purpose, we introduce K — a set of clusters k to group observations x into. We also assume each cluster has an unknown prior probability p(k) and that observations and clusters are statistically connected by an unknown joint probability p(x, k). Let us apply the conditional probability formula in reverse to expand pᵗ(x):
Then, for an arbitrary k such that pᵗ(k | x) and pᵗ(k) are non-zero, we have
We have a problem: if the posterior probability in the denominator is zero, the entire sum is undefined (division by zero). For each x, we need to use only such k that avoids this uncanny situation. I will introduce predicates, π, that say whether we can use the specific term or whether it yields an incorrect expression, and define that, if π is zero, the logarithm is ignored:
How do we unambiguously choose which predicate is non-zero? Before answering this question, I will relax the conditions and use "fuzzy" predicates — weights:
Now, let us think of a non-negative quantity that depends on k and x, sums to 1 over all elements of K, and becomes 0 whenever the joint probability of x and k is zero. It is the posterior probability!
Can the numerator be zero when the posterior probability is not? No. Because the probability of any observation is non-zero (if we faced the observation x, it is probable, thus p(x) > 0), and the implication is just an application of the conditional probability formula:
By applying the conditional probability formula once more, we obtain the numerator, proving it is greater than zero if the posterior probability is greater than zero:
Improving The Log-Likelihood
We aim to obtain probabilities that increase the log-likelihood, which is the iteration of the EM algorithm. I will expand the logarithm to reveal that the log-likelihood contains an entropy and cross-entropies:
Prior Probabilities
As long as prior probabilities pᵗ(k) are not tied to any statistical model, we can optimise the member that contains them without significant complications. To solve the problem, we can use the method of Lagrange multipliers, but I will go straight to the solution:
Conditional Probabilities
Optimising the first sum in the expanded log-likelihood expression is less trivial to achieve. If we used the method of Lagrange multipliers again, we would achieve nothing good, because by Gibbs' inequality, the answer would be p(x | k) = p(k | x). We must pick constraints for the probabilities — a [generative] statistical model P, and pick such a probability from it that does not decrease the sum:
When possible, we want to solve the optimisation problem
Posterior Probabilities
The last sum in the expanded log-likelihood expression contains logarithms of conditional probabilities that can be computed using Bayes' rule from the probabilities we already altered, so we apply the rule:
The Gibbs' inequality states that the information entropy of a distribution is not greater than its cross-entropy with any other distribution. The modification of posterior probabilities changed the entropy into a cross-entropy, so it is not less than the original value:
Monotonicity of The EM Algorithm
Each summand of the expanded log-likelihood did not decrease after the algorithm iteration, which means the entire expression did not decrease:
Correctness of The EM Algorithm
Using the conditional probability formulae as we did in the beginning, we fold the right-hand side of the inequality:
The logarithm does not depend on k, so it can be moved out of the inner sum using the distributivity of multiplication over addition:
Then the inner sum equals 1, because it sums probabilities over the entire set of independent events, so we get our new log-likelihood
We have what we wanted — the quantity that does not decrease after any algorithm iteration is the log-likelihood of the input sample.
The EM Algorithm
Input
- A sample, x, of observations (order and repetitions do not matter).
- A finite number, |K| > 2, of clusters.
- A non-empty family, P ∋ pₓ|ₖ, of generative models.
Initialisation
I described “the next” step of the algorithm, but I never mentioned the first step. Honestly, I do not know any general good advice on initialisation, and I guess nobody does (maybe a general simple recommendation does not exist). What I can advise is: if you do not have anything better, initialise the clusters randomly.
Get Valeriy Krygin’s stories in your inbox
Join Medium for free to get updates from this writer.
For example, for each x, you can select a random k, and say p⁰(k | x) = 1, and, for all other k’ ≠ k, say p⁰(k’ | x) = 0. From here, you can compute p¹(k) and p¹(x | k), and proceed. This may work; sometimes it may not.
Alternatively, you may sample parameters of the generative model, p⁰(x | k), and the prior probabilities, p⁰(k), but the exact initial values sampling method (or hand picking) depends on your problem.
Expectation (The E-Step)
This is the easy part of the algorithm, but computationally nasty because it involves averaging of potentially very small numbers:
Maximisation (The M-Step)
In general, this is a very hard part, because maximisation problems can be very hard:
Result
The result is an improved (or at least not degraded) log-likelihood of the input sample:
The estimated probabilities constitute the interesting part of the result: we can use them for visualisation and even future classification.
Conclusions and Further Directions
I mentioned a generative statistical model, P. The algorithm does not specify its nature, so the Gaussian Mixture Model is not the only model it applies to. P can be a family of probabilistic graphical models — then the maximisation becomes intractable, but small improvements are computationally feasible.
Here is the homework for you:
- Prove that applying Bayes' rule does not introduce infinities into the expression.
- Implement the EM Algorithm for the classical Gaussian Mixture Model example. What happens if multiple clusters share the same centroid?
- Implement a generic EM Algorithm that accepts an arbitrary generative model and maximisation procedure. How do you handle the machine-zero issue when averaging generally?
- Apply the EM Algorithm to a problem involving Markov Chains (so the Maximisation step requires dynamic programming) — this particular case of the algorithm is also known as the Baum-Welch algorithm.
- Read Ten Lectures on Statistical and Structural Pattern Recognition (or at least The interaction of learning and self-organization in pattern recognition).
- Find and understand the proofs of the method of Lagrange multipliers and Gibbs’ inequality.
Good luck, and have fun!
