Due: Friday, Nov. 13.
Every estimate you have computed in this course used data you could see. The EM algorithm computes maximum likelihood estimates when part of the data is missing, and here the missing part is the answer key. Each data point came from one of several Gaussian clusters, and nobody records which. EM estimates the cluster parameters anyway, which makes it a clustering algorithm. Then you will use it to answer a harder question that has no obvious likelihood at all: how many clusters are there?
The lab has two parts. In Part 1 you fit the clusters when you already know how many there are. In Part 2 you let the data tell you how many clusters there are. Both parts turn out to be the same procedure, seen from two sides, and finding that out is the point of the lab.
This lab also teaches something about building software with an AI. In the earlier labs the page handed you the exact signature of every function. This one does not, and that is deliberate. Here the program has real structure, so you will design it yourself: for each piece you write the specification, saying what it computes, what goes in, and what comes out, and only then does the AI write the code. Writing a specification precise enough that a stranger could implement it is the skill this lab grades.
The background is Chapter 12. The notation here matches the book.
labs conda environment from Lab 2. All code
in python, using numpy or PyTorch.rng = np.random.default_rng(0). State in
your report which generator you used.image-processing-labs
repository, in a directory named lab6.Let \( \{X_n\}_{n=1}^{N} \) be i.i.d. discrete labels taking values in \( \{0, \cdots, M-1\} \) with \( P\{X_n = m\} = \pi_m \). Let \( \{Y_n\}_{n=1}^{N} \), with \( Y_n \in \mathbb{R}^p \), be conditionally independent Gaussian random vectors given the labels, with \( Y_n \sim N(\mu_{x_n}, R_{x_n}) \). Write \( \theta_m = (\mu_m, R_m) \) for the parameters of component \( m \), so that
$$ p(y_n | \theta_m) = \frac{1}{(2\pi)^{p/2} |R_m|^{1/2}} \exp\left\{ -\frac{1}{2} (y_n - \mu_m)^t R_m^{-1} (y_n - \mu_m) \right\} , \tag{1} $$and the full parameter vector is \( \theta = (\pi_0, \mu_0, R_0, \cdots, \pi_{M-1}, \mu_{M-1}, R_{M-1}) \). In your code, keep \( \theta \) as three arrays: the weights of shape \( (M,) \), the means of shape \( (M, p) \), and the covariances of shape \( (M, p, p) \).
The labels are never observed. Summing them out gives the Gaussian mixture distribution
$$ p(y_n | \theta) = \sum_{m=0}^{M-1} \pi_m \, p(y_n | \theta_m) , \tag{2} $$and, since the samples are i.i.d., the log likelihood of the whole data set is
$$ \log p(y | M, \theta) = \sum_{n=1}^{N} \log \left( \sum_{m=0}^{M-1} \pi_m \, p(y_n | \theta_m) \right) . \tag{3} $$We want the \( \theta \) that maximizes (3). Direct maximization is hard, and the reason is visible in the formula: there is a sum inside the logarithm, so the log does not distribute and nothing separates.
Suppose for a moment that the labels \( x_n \) were observed. Then estimation would be easy, and it is worth seeing how easy, because EM is built directly on this case.
Everything the ML estimate needs to know about labeled data is contained in three statistics per class:
$$ N_m = \sum_{n=1}^{N} \delta(x_n - m) , \qquad b_m = \sum_{n=1}^{N} y_n \, \delta(x_n - m) , \qquad S_m = \sum_{n=1}^{N} y_n y_n^t \, \delta(x_n - m) , \tag{4} $$where \( \delta(\cdot) \) is 1 at zero and 0 elsewhere. So \( N_m \) counts the points in class \( m \), \( b_m \) sums them, and \( S_m \) sums their outer products. These three are called the sufficient statistics of the labeled data, because the ML estimate depends on the data only through them. Chapter 12 shows that the ML estimates are simple functions of them:
$$ \hat{\pi}_m = \frac{N_m}{N} , \qquad \hat{\mu}_m = \frac{b_m}{N_m} , \qquad \hat{R}_m = \frac{S_m}{N_m} - \frac{b_m b_m^t}{N_m^2} . \tag{5} $$These are just the class fraction, the class sample mean, and the class sample covariance. Keep (4) and (5) in view for the next step.
The labels are not observed, so the statistics (4) cannot be computed. The EM idea is one sentence: replace each statistic by its expected value given the observed data and the current parameter estimate, then apply (5) as if the expected statistics were the real thing.
The expectation of \( \delta(X_n - m) \) given \( Y = y \) is the posterior probability that sample \( n \) belongs to class \( m \). The E-step computes it for every \( n \) and \( m \):
$$ P_{n,m} \leftarrow P\{ X_n = m \,|\, Y = y, \theta \} = \frac{ \pi_m \, p(y_n | \theta_m) } { \displaystyle \sum_{j=0}^{M-1} \pi_j \, p(y_n | \theta_j) } . \tag{6} $$These are the soft assignments. Sample \( n \) has partial membership \( P_{n,m} \) in cluster \( m \), and \( \sum_m P_{n,m} = 1 \). The expected statistics are then (4) with \( \delta(x_n - m) \) replaced by \( P_{n,m} \):
$$ \bar{N}_m \leftarrow \sum_{n=1}^{N} P_{n,m} , \qquad \bar{b}_m \leftarrow \sum_{n=1}^{N} y_n P_{n,m} , \qquad \bar{S}_m \leftarrow \sum_{n=1}^{N} y_n y_n^t P_{n,m} , \tag{7} $$and the M-step is (5) applied to the expected statistics:
$$ \pi_m \leftarrow \frac{\bar{N}_m}{N} , \qquad \mu_m \leftarrow \frac{\bar{b}_m}{\bar{N}_m} , \qquad R_m \leftarrow \frac{\bar{S}_m}{\bar{N}_m} - \frac{\bar{b}_m \bar{b}_m^t}{\bar{N}_m^2} . \tag{8} $$One EM iteration is one E-step followed by one M-step. Chapter 12 proves the property that makes the algorithm usable: each iteration increases the log likelihood (3), or leaves it unchanged. It never decreases. That property is also the best test you have of whether your code is right, and you will check it numerically in Step 7.
Generate \( N = 500 \) samples from a Gaussian mixture with \( p = 2 \), \( M = 3 \), and \( \pi = [0.4, 0.4, 0.2] \), means
$$ \mu_0 = \begin{bmatrix} 2 \\ 2 \end{bmatrix} , \qquad \mu_1 = \begin{bmatrix} -2 \\ -2 \end{bmatrix} , \qquad \mu_2 = \begin{bmatrix} 5.5 \\ 2 \end{bmatrix} , \tag{9} $$and covariances
$$ R_0 = \begin{bmatrix} 1 & 0.1 \\ 0.1 & 1 \end{bmatrix} , \qquad R_1 = \begin{bmatrix} 1 & -0.1 \\ -0.1 & 1 \end{bmatrix} , \qquad R_2 = \begin{bmatrix} 1 & 0.2 \\ 0.2 & 0.5 \end{bmatrix} . \tag{10} $$
Use rng = np.random.default_rng(0). Draw each label
\( x_n \) from \( \pi \), then draw \( y_n \) from
\( N(\mu_{x_n}, R_{x_n}) \). Save both \( y \), of shape
\( 500 \times 2 \), and the true labels \( x \) to a file, in the
order they were generated. You will reuse this exact data set in every
remaining step, so generate it once and load it after that.
The parameters above fully define the distribution you are sampling from. In Parts 1 and 2 you will try to recover this same \( \theta \) from the samples alone, without the labels.
A scatter plot of the 500 samples, with the three true classes in three colors and a legend. Use equal axis scaling, so that the covariances are not distorted by the aspect ratio.
Do this step without AI help. Both derivations are short, and they are what make the rest of the lab make sense rather than being a sequence of formulas to type in.
Both derivations, handwritten or typeset, a few lines each.
Now design the implementation. Do not ask the AI for "an EM program." Build it from the four functions whose signatures are given below. The parameters are the three arrays from Step 1: the weights, of shape \( (M,) \), the means, of shape \( (M, p) \), and the covariances, of shape \( (M, p, p) \). Every function reads or returns them. Assembled, the four functions are a complete EM program.
The signatures are fixed, so that everyone's code has the same structure. Before you write any code for a function, write its specification, which has three parts:
def e_step(Y, weights, means, covariances):
# posterior class probabilities P_{n,m}, equation (6)
# Y (N,p), weights (M,), means (M,p), covariances (M,p,p)
# returns P, shape (N,M)
def m_step(Y, P):
# expected statistics (7) and updated parameters (8)
# Y (N,p), P (N,M); add 1e-6 * I to each covariance
# returns weights (M,), means (M,p), covariances (M,p,p)
def log_likelihood(Y, weights, means, covariances):
# the log likelihood of equation (3); returns a float
def run_em(Y, weights, means, covariances, num_iterations):
# num_iterations of e_step then m_step
# returns weights, means, covariances, and the list of log
# likelihoods recorded after each iteration
Here is how the four functions fit together in one EM run. This is the
body of run_em, and it is the whole algorithm.
Given data Y and initial weights, means, covariances:
record log_likelihood(Y, weights, means, covariances)
repeat 20 times:
P = e_step(Y, weights, means, covariances)
weights, means, covariances = m_step(Y, P)
record log_likelihood(Y, weights, means, covariances)
The four specifications, verbatim. These are a graded part of the lab, not a preliminary.
Hand the AI one specification at a time, and test each function before you build the next.
Two implementation requirements apply to everyone, and both exist because of failures that produce no error message.
Run these four checks.
e_step. Every row of \( P_{n,m} \) sums to
one, to within numerical precision. Include a test point placed
far from all component means, where an implementation that skipped
the log domain returns NaN.m_step. Feed it the hard assignments
\( P_{n,m} = \delta(x_n - m) \) built from the true labels of
Step 4. By your second derivation in Step 5, the output must equal
the class fractions, class sample means, and class sample
covariances computed directly from the labeled data. Verify that
it does.log_likelihood. For \( M = 1 \), the mixture
(2) is a single Gaussian. Check the function against a direct
evaluation of the Gaussian log density summed over the
samples.run_em. The recorded log likelihood is
nondecreasing from one iteration to the next. Report the largest
decrease you observe, which should be zero or a rounding-level
negative number. Also check that if you initialize at the true
parameters of Step 4, the parameters stay close to the true
values.The result of each of the four checks, with the numbers.
Run EM on the data from Step 4, with \( M = 3 \), for 20 iterations, and the initialization
$$ \pi \leftarrow [1/3, 1/3, 1/3] , \qquad R_m \leftarrow \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \ \mbox{ for } m = 0, 1, 2 , \tag{11} $$with \( \mu_0, \mu_1, \mu_2 \) set to the first three sample vectors produced in Step 4.
We will track the cost, which is the negative log likelihood \( -\log p(y | M, \theta) \). EM increases the log likelihood, so it decreases the cost. Plot the cost against the iteration number and look at the last few iterations: the cost is still dropping at iteration 20, so EM has not finished. The fit is still usable. With the right number of clusters, 20 iterations buys a reasonable answer, not a great one. More iterations would keep improving it.
A plot of the cost \( -\log p(y | M, \theta) \) versus iteration number, for iterations 1 through 20. State whether the cost is still decreasing at iteration 20.
A table of the estimated \( \pi_m \), \( \mu_m \), and \( R_m \), each placed next to the true value it corresponds to, using the correspondence you found.
A scatter plot of the data colored by the hard label \( \hat{x}_n = \arg\max_m P_{n,m} \), with the estimated mean of each component marked, and the fraction of samples whose hard label disagrees with the true label from Step 4. The hard label is the most-probable component, the one with the largest posterior \( P_{n,m} \). It uses the weights \( \pi_m \), so it is not in general the component with the largest density.
In Step 8 you set the initial means to the first three samples. That is one choice. The log likelihood of a mixture has many local maxima, and EM goes uphill to the nearest one, so a different start can reach a different answer.
Repeat Step 8 five times with \( M = 3 \), keeping \( \pi \) and \( R_m \) as in (11), but choosing the three initial means as three distinct data points drawn at random from the 500 samples. For run \( s \), use
idx = np.random.default_rng(s).choice(500, size=3, replace=False) initial_means = y[idx]
for \( s = 0, 1, 2, 3, 4 \). Run 20 iterations each.
One plot with the five cost curves plus the curve from Step 8, on the same axes with a legend, and a table of the final cost for all six runs. State whether your Step 9 prediction was right.
Five scatter plots, one per random start, each colored by the hard label and marked with the estimated means. For each, give the final cost. In one or two sentences, say which starts reached a good clustering and which did not, and what the poor ones did with the clusters.
So far you told the algorithm that \( M = 3 \). In practice nobody tells you. And you cannot choose \( M \) by maximizing the likelihood, because any order-\( M \) model is a special case of an order-\( (M+1) \) model, so the maximized likelihood can only go up with \( M \). Left alone, the likelihood will happily give every data point its own cluster.
The fix is to add a penalty that grows with the number of parameters. The minimum description length criterion is
$$ MDL(M, \theta) = -\log p(y | M, \theta) + \frac{1}{2} L(M) \log(Np) , \tag{12} $$where \( L(M) \) is the number of real numbers needed to specify \( \theta \). Each component needs 1 number for \( \pi_m \), \( p \) numbers for \( \mu_m \), and \( (p+1)p/2 \) numbers for the symmetric matrix \( R_m \), and one number is redundant because the \( \pi_m \) sum to one. So
$$ L(M) = M \left( 1 + p + \frac{(p+1)p}{2} \right) - 1 . \tag{13} $$The first term of (12) is exactly the cost you plotted in Part 1. So for a fixed number of clusters \( M \), the penalty \( \tfrac12 L(M) \log(Np) \) is a constant: it does not depend on \( \theta \). The MDL is then the cost plus that constant, and minimizing MDL over \( \theta \) at fixed \( M \) is the same as maximizing the likelihood, which is exactly what EM already does. The penalty changes nothing within a single \( M \). Its only job is to compare different values of \( M \): more clusters means more parameters, means a larger constant, and that is what stops the search from giving every point its own cluster.
To find \( \hat{M} \), start with more components than you expect and remove them one at a time.
Order estimation by merging
Initialize M_o to the largest order you will consider, and
initialize theta
For M = M_o down to 1:
Run EM to produce theta*(M)
Record MDL(M) using (12)
Merge two components of theta*(M) to initialize order M-1
Select M_hat as the value of M minimizing MDL(M)
Merge the pair of components \( (l, m) \) that minimizes the distance
$$ d(l,m) = \frac{N \pi_l}{2} \log\left( \frac{|R_{(l,m)}|}{|R_l|} \right) + \frac{N \pi_m}{2} \log\left( \frac{|R_{(l,m)}|}{|R_m|} \right) , \tag{14} $$where the merged component has parameters
$$ \pi_{(l,m)} = \pi_l + \pi_m , \qquad \mu_{(l,m)} = \frac{\pi_l \mu_l + \pi_m \mu_m}{\pi_l + \pi_m} , \tag{15} $$ $$ \begin{aligned} R_{(l,m)} = \frac{1}{\pi_l + \pi_m} \Big[ &\pi_l \left( R_l + (\mu_l - \mu_{(l,m)})(\mu_l - \mu_{(l,m)})^t \right) \\ &+ \pi_m \left( R_m + (\mu_m - \mu_{(l,m)})(\mu_m - \mu_{(l,m)})^t \right) \Big] . \end{aligned} \tag{16} $$This distance upper bounds the increase in MDL caused by the merge, so merging the closest pair gives up as little as possible at each reduction in order. After the merge, the merged parameters become the new component and the others carry over unchanged.
Each lower order starts from the merge of the order above it. This is a warm start: the order above has already spread its components across the data, so after a merge the remaining components sit on real cluster regions, and EM finishes in a few iterations. That is why 20 iterations per order is more than enough here.
Three more functions are needed, with the signatures below. Specify each one the same way you did in Step 6, before any code.
def mdl(Y, weights, means, covariances):
# the MDL value of equations (12) and (13); returns a float
def merge_closest(Y, weights, means, covariances):
# finds the closest pair (14), merges it with (15) and (16)
# returns the order M-1 weights, means, covariances
def estimate_order(Y, max_order, initial_means, num_iterations):
# the search loop: for M from max_order down to 1, run_em then
# merge_closest, recording mdl(M) at each order
# returns, for each M, the MDL value and the parameters, and the
# per-iteration MDL values for the D11 plot
Run it on the data from Step 4 with \( M_o = 9 \), the initialization \( \pi \leftarrow [1/9, \cdots, 1/9] \), \( R_m \leftarrow I \) for all \( m \), and \( \mu_0 \) through \( \mu_8 \) set to the first nine sample vectors from Step 4. Run 20 EM iterations for each value of \( M \).
The three additional specifications, verbatim.
A plot of the MDL value versus the cumulative number of EM iterations over the whole run. Plot the MDL value after each EM iteration, and when an iteration is followed by a merge, also plot the value immediately after the merge, so that merge steps show two points. Mark the merge locations.
A plot of \( MDL(M) \) versus \( M \) for \( M = 1, \cdots, 9 \), with the minimizing \( \hat{M} \) labeled, and a statement of whether \( \hat{M} \) equals the true value \( M = 3 \).
The merge search passed through an \( M = 3 \) model on its way down. Compare its cost there with the five \( M = 3 \) costs from Step 9. The merge search reaches a good answer, the same quality as the best random start, and it does so without being told the number of clusters. Starting with too many components and merging down is not only how you count the clusters. It is also a reliable way to fit them, because it reaches the good answer that a single random start often misses.
Answer each question in a few sentences, using your own plots and numbers as evidence.
Your answers to the six questions above.
Submit two files through Brightspace.
.py files. Do not paste the code
into the report and do not send a zip file.
Also commit your code to your image-processing-labs
repository under lab6 and push it before you submit.
Concise and clear beats long.
EM is the standard tool for fitting models with hidden variables, and the pattern you used here, expected sufficient statistics in place of observed ones, carries far beyond Gaussian mixtures. The deeper result of this lab is that fitting the clusters, avoiding a bad answer from a poor starting point, and choosing the number of clusters are one procedure: start with too many components and merge down. You met it as three separate problems and found it was one.