Gaussian Mixtures and the EM Algorithm

A Gaussian mixture assumes each point was generated by one of \(K\) Gaussian components, and that the component label was not observed. This missing label is a latent variable. Expectation–maximization (EM) fits the model by alternating between estimating those labels probabilistically and updating the component parameters. EM also applies to other latent-variable models, including hidden Markov models; some models require approximate versions of its steps.

The model and the two steps

The density is \(p(x)=\sum_k\pi_k\mathcal N(x\mid\mu_k,\Sigma_k)\), where \(\pi_k\geq0\), \(\sum_k\pi_k=1\), and each covariance \(\Sigma_k\) is positive definite. The weight \(\pi_k\) is the probability of selecting component \(k\) before drawing a point from it. With observed component labels and separate covariances, we could estimate each Gaussian from its assigned points and each weight from its fraction of the sample. EM uses fractional assignments because those labels are missing.

• E-step. Compute each point’s responsibility \(r_{ik} = \pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k) / \sum_j \pi_j \mathcal{N}(x_i \mid \mu_j, \Sigma_j)\), conditional on the current fitted parameters.
• M-step. Let \(N_k=\sum_i r_{ik}\), the component’s fractional sample count. Update \(\pi_k=N_k/n\), \(\mu_k=\sum_i r_{ik}x_i/N_k\), and \(\Sigma_k=\sum_i r_{ik}(x_i-\mu_k)(x_i-\mu_k)^\top/N_k\), using the newly updated mean. These are maximum-likelihood updates, so the covariance denominator is \(N_k\).

For example, suppose two weights are 0.6 and 0.4, and their densities at one point are 0.2 and 0.1. Their contributions are 0.12 and 0.04, so the responsibilities are \(0.12/0.16=0.75\) and 0.25. That point contributes three quarters of an observation to the first component’s update. In the code, R has one row per observation and one column per component; every row sums to one.

The total log-likelihood is \(\ell=\sum_i\log p(x_i)\), which adds the log density of every observation. Larger values indicate a better density fit on that same sample. score(X) returns its per-observation average, so the code multiplies by n. Natural-log units are called nats. Run the body blocks in order because later comparisons reuse the first dataset X.

import numpy as np
from scipy.stats import multivariate_normal
from scipy.special import logsumexp
from sklearn.cluster import kmeans_plusplus
from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(0)
true_mu = np.array([[0, 0], [4, 1], [1, 5]], float)
true_cov = [np.array([[1, .6], [.6, 1]]), np.array([[1.5, -.5], [-.5, .7]]),
            np.array([[.6, 0], [0, 2.0]])]
n, K = 3000, 3
comp = rng.choice(3, n, p=[.40, .35, .25])
X = np.vstack([rng.multivariate_normal(true_mu[k], true_cov[k]) for k in comp])

mu, _ = kmeans_plusplus(X, K, random_state=0)
cov = [np.cov(X.T) for _ in range(K)]
pi = np.full(K, 1 / K)
initial_mu, initial_cov, initial_pi = mu.copy(), np.array(cov), pi.copy()
history = []
for it in range(1001):
    logp = np.column_stack([np.log(pi[k])
                           + multivariate_normal(mu[k], cov[k]).logpdf(X)
                           for k in range(K)])
    log_density = logsumexp(logp, axis=1)
    history.append(log_density.sum())
    if it:
        gain = (history[-1] - history[-2]) / n
        if gain < -1e-12:
            raise RuntimeError("Likelihood decreased beyond roundoff tolerance")
        if abs(gain) < 1e-10:
            break
    if it == 1000:
        raise RuntimeError("Iteration limit reached")
    R = np.exp(logp - log_density[:, None])
    Nk = R.sum(axis=0)
    if np.any(Nk <= 0):
        raise RuntimeError("An empty component needs reinitialization")
    pi = Nk / n
    mu = (R.T @ X) / Nk[:, None]
    cov = [(R[:, k:k+1] * (X - mu[k])).T @ (X - mu[k]) / Nk[k]
           for k in range(K)]

sk = GaussianMixture(3, n_init=10, init_params="k-means++", random_state=0).fit(X)
matched = GaussianMixture(3, n_init=1, weights_init=initial_pi,
    means_init=initial_mu, precisions_init=np.linalg.inv(initial_cov),
    reg_covar=0, tol=1e-10, max_iter=1000, random_state=0).fit(X)
print(f"scratch M-steps {it}, log-likelihood {history[-1]:.6f}")
print(f"sklearn default settings          {sk.score(X) * n:.6f}")
print(f"sklearn matched start and tol     {matched.score(X) * n:.6f}")
print(f"monotone within tolerance: {all(np.diff(history) >= -1e-9)}")
print(f"first five: {np.round(history[:5], 2)}")
# scratch M-steps 37, log-likelihood -11131.124979
# sklearn default settings          -11132.245573
# sklearn matched start and tol     -11131.124979
# monotone within tolerance: True
# first five: [-13626.02 -11800.62 -11297.9  -11162.69 -11141.55]

The scratch implementation reaches \(-11131.124979\) after 37 M-steps. The library fit with default tolerance stops at \(-11132.245573\). That difference alone does not establish different local optima: the fits also use different stopping rules, initialization, and covariance regularization. Matching the initial parameters, setting regularization to zero, and tightening the library tolerance gives the same printed likelihood as the scratch implementation.

With an exact E-step and an M-step that maximizes its expected complete-data log-likelihood, EM cannot decrease the observed-data likelihood in exact arithmetic. The E-step constructs a lower bound that touches the current likelihood, and the M-step raises that bound. A material decrease in this unregularized implementation needs investigation; tiny floating-point changes, approximate steps, and changes to the objective require separate interpretation. Nondecreasing likelihood does not guarantee a global maximum, or even that the parameters approach a finite local maximum.

The E-step uses logsumexp to evaluate the logarithm of a sum without first forming very small densities. It effectively subtracts the largest log value before exponentiating and adds it back afterward. Direct density calculations can underflow in high dimensions or far into a distribution’s tails; there is no fixed dimension at which this starts. This teaching implementation leaves covariance regularization out so the likelihood guarantee applies directly. It assumes nonsingular updates; the collapse example below explains why production code needs additional safeguards.

How k-means relates to a mixture

With shared covariance \(\Sigma_k=\sigma^2 I\), fixed positive mixing weights, and a uniquely nearest mean, the responsibility of that mean approaches one as \(\sigma\to0\). A tie between nearest means retains their relative weights. Hard assignments followed by averages of the assigned points recover Lloyd’s updates, provided empty clusters are handled.

Another connection is exact at fixed variance: hard-assignment likelihood with equal mixture weights and a fixed shared spherical covariance is maximized by minimizing the k-means sum of squared distances. Finite-variance mixture EM still uses soft assignments. In the library, spherical allows a different spherical variance for each component; tied shares one full covariance. The diag option gives each component axis-aligned ellipsoids, while full permits its own tilted ellipsoid. These covariance families do not form a single nested sequence.

import numpy as np
from sklearn.mixture import GaussianMixture

for ct in ("spherical", "diag", "tied", "full"):
    g = GaussianMixture(3, covariance_type=ct, n_init=10, tol=1e-7, max_iter=1000,
                        init_params="k-means++", random_state=0).fit(X)
    d = X.shape[1]
    covariance_params = {"spherical": 3, "diag": 3*d,
                         "tied": d*(d+1)//2, "full": 3*d*(d+1)//2}
    params = 3*d + 2 + covariance_params[ct]
    print(f"{ct:10s} params {params:4d}"
          f"  log-likelihood {g.score(X) * n:10.1f}  BIC {g.bic(X):10.1f}")
g_tied = GaussianMixture(3, covariance_type="tied", n_init=10,
    init_params="kmeans", tol=1e-7, max_iter=1000, random_state=0).fit(X)
print(f"tied with kmeans initialization: log-likelihood {g_tied.score(X)*n:.1f},"
      f" BIC {g_tied.bic(X):.1f}")
# spherical  params   11  log-likelihood   -11731.6  BIC    23551.3
# diag       params   14  log-likelihood   -11410.8  BIC    22933.7
# tied       params   11  log-likelihood   -13194.9  BIC    26477.8
# full       params   17  log-likelihood   -11131.1  BIC    22398.4
# tied with kmeans initialization: log-likelihood -11730.5, BIC 23549.0

The full fit has the largest likelihood and smallest BIC among these fits. Its covariance family includes the data-generating shapes, but a fitted score also depends on optimization. In particular, the poor tied result with k-means++ initialization should not be attributed entirely to sharing the covariance. The additional kmeans initialization check obtains a much larger likelihood within the same tied-covariance family.

For \(K\) components in \(d\) dimensions, the covariance parameter counts are \(K\) for spherical, \(Kd\) for diagonal, \(d(d+1)/2\) for tied, and \(Kd(d+1)/2\) for full. Add \(Kd\) mean parameters and \(K-1\) free mixing weights to obtain the totals printed above. At \(d=50,K=5\), the covariance counts alone are 5, 250, 1275, and 6375. More parameters require stronger data support, but soft assignments do not give a simple integer “rows per component” rank rule. A weighted covariance can have full rank when the centered, positively weighted points span the space, yet be unstable when most weights are tiny. Inspect component mass, covariance eigenvalues, and sensitivity to regularization as well as the parameter count.

The likelihood is unbounded

For an unconstrained mixture with at least two components and separate covariances, the sample likelihood is unbounded. Give one component positive weight, center it on a data point, and shrink its covariance toward zero. Its density at that point diverges while another component can retain positive density at every remaining point. Thus this likelihood has no finite global maximizer. Covariance restrictions, regularization, or suitable priors change the fitting problem.

import numpy as np
from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(0)
Xd = np.vstack([rng.normal([0, 0], 1.0, (200, 2)),
                np.tile([3.0, 3.0], (4, 1))])   # four identical points

for reg in (1e-10, 1e-6):
    best = None
    for seed in range(60):
        g = GaussianMixture(6, covariance_type="full", reg_covar=reg, n_init=1,
                            init_params="random", random_state=seed,
                            max_iter=1000).fit(Xd)
        ll = g.score(Xd) * len(Xd)
        if best is None or ll > best[0]:
            best = (ll, g)
    ll, g = best
    smallest = min(np.linalg.det(c) for c in g.covariances_)
    print(f"reg_covar {reg:.0e}  max log-likelihood {ll:10.2f}"
          f"   smallest component determinant {smallest:.3e}")
# reg_covar 1e-10  max log-likelihood    -491.74   smallest component determinant 1.000e-20
# reg_covar 1e-06  max log-likelihood    -528.58   smallest component determinant 1.000e-12

The smallest determinant is approximately \(\text{reg\_covar}^2\): \(10^{-20}\) and \(10^{-12}\) in two dimensions. In scikit-learn, reg_covar is added to each covariance diagonal after the weighted covariance calculation. Here a component concentrates on the repeated point, making its raw covariance nearly zero, so the stored covariance is approximately \(\text{reg\_covar}I\). Its size is set by the covariance update; the separate likelihood-change tolerance controls when iteration stops.

Reducing regularization raises the best fitted log-likelihood from \(-528.58\) to \(-491.74\). Together with the duplicated points and near-zero raw covariance, this is evidence of concentration on those points. An increase in likelihood alone would not diagnose collapse. A fit can meet its stopping criterion for each fixed regularization value even though allowing that value to approach zero exposes the unbounded likelihood.

Inspect component weights and covariance eigenvalues alongside the observations receiving high responsibility. Determinants depend strongly on dimension and feature units, so a small determinant alone does not diagnose collapse. Increasing reg_covar, reducing \(K\), or restricting covariances are options to compare using the intended evaluation criterion. Covariance priors can discourage collapse, with behavior depending on the prior and inference method. The GaussianMixture documentation describes the library’s regularization and stopping parameters.

Choosing K

import numpy as np
from sklearn.mixture import GaussianMixture

print(f"{'K':>3} {'BIC':>12} {'AIC':>12}")
for K_ in range(1, 8):
    g = GaussianMixture(K_, covariance_type="full", n_init=10, tol=1e-7, max_iter=1000,
                        init_params="k-means++", random_state=0).fit(X)
    print(f"{K_:3d} {g.bic(X):12.1f} {g.aic(X):12.1f}")
#   K          BIC          AIC
#   1      26429.7      26399.7
#   2      23591.3      23525.2
#   3      22398.4      22296.3
#   4      22430.5      22292.4
#   5      22465.4      22291.2
#   6      22506.9      22296.7
#   7      22545.2      22299.0

For these fitted candidates, BIC selects \(K=3\), matching the generating component count. AIC selects \(K=5\), with 22291.2 against 22296.3 for three components. Its lighter penalty favors extra components in this sample. These criteria compare the optimized fits obtained here, so inadequate optimization can also affect the ranking.

Writing \(\ell\) for fitted total log-likelihood and \(p\) for the number of free parameters, \(\mathrm{AIC}=-2\ell+2p\) and \(\mathrm{BIC}=-2\ell+p\log n\); smaller values are preferred. With \(n=3000\), BIC’s penalty is about \(8p\), four times AIC’s. AIC is motivated by predictive fit, while BIC is often used for component selection. Finite mixtures require additional care with these asymptotic arguments: redundant components and collapsing covariances violate assumptions used in the usual BIC derivation. Drton and Plummer discuss these singular-model complications. Consistency claims for selecting a true component count need additional identifiability and parameter-space assumptions; merely including the true model among candidates is insufficient.

When no finite Gaussian mixture adequately describes the population, extra components can approximate shape details without identifying additional meaningful groups. Even a good density fit need not have one component per substantive group. Use the BIC comparison together with fit diagnostics and the purpose of the model; it does not establish how many real-world categories exist.

Exercises

1. Soft assignments and abstention. Fit two overlapping components on training data. On independent test data, measure the fraction retained and the actual label accuracy at several responsibility thresholds. Also report ARI, keeping it distinct from accuracy.

The synthetic component labels let you evaluate the fit. Match component IDs to labels using training data only, then keep that mapping fixed for test evaluation.

import numpy as np
from scipy.optimize import linear_sum_assignment
from sklearn.mixture import GaussianMixture
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(1)
def sample(n):
    y = rng.integers(0, 2, size=n)
    return rng.normal(size=(n, 2)) + np.column_stack((2.2*y, np.zeros(n))), y
X_tr, y_tr = sample(1000)
X_te, y_te = sample(4000)
g = GaussianMixture(2, n_init=10, init_params="k-means++",
                    tol=1e-7, max_iter=1000, random_state=0).fit(X_tr)
km = KMeans(2, n_init=10, random_state=0).fit(X_tr)
train_labels = g.predict(X_tr)
counts = np.zeros((2, 2), dtype=int)
np.add.at(counts, (train_labels, y_tr), 1)
rows, cols = linear_sum_assignment(-counts)
mapping = np.empty(2, dtype=int)
mapping[rows] = cols
labels = g.predict(X_te)
pred = mapping[labels]
confidence = g.predict_proba(X_te).max(axis=1)
print(f"test GMM ARI {adjusted_rand_score(y_te, labels):.4f}"
      f"   k-means ARI {adjusted_rand_score(y_te, km.predict(X_te)):.4f}")
for t in (0.50, 0.70, 0.90, 0.99):
    keep = confidence >= t
    accuracy = (pred[keep] == y_te[keep]).mean()
    print(f"threshold {t:.2f}: coverage {keep.mean():6.1%}, accuracy {accuracy:.4f}")
# test GMM ARI 0.5572   k-means ARI 0.5624
# threshold 0.50: coverage 100.0%, accuracy 0.8732
# threshold 0.70: coverage  84.9%, accuracy 0.9161
# threshold 0.90: coverage  61.2%, accuracy 0.9653
# threshold 0.99: coverage  23.6%, accuracy 0.9936

2. Mixtures for density estimation. Fit mixtures to a uniform square, which a finite nonsingular Gaussian mixture cannot represent exactly. Compare held-out log densities with a Gaussian kernel density estimate.

Compare the average log density near the edge with the interior. Read the reported values as log densities and compare them with the known uniform density.

import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import KernelDensity

rng = np.random.default_rng(0)
X_tr = rng.uniform(-1, 1, (4000, 2))           # a uniform square
X_te = rng.uniform(-1, 1, (4000, 2))
truth = np.log(0.25)                           # density of U([-1,1]^2) is 1/4

print(f"true log-density {truth:.4f}")
for K in (1, 4, 16, 64):
    g = GaussianMixture(K, n_init=3, init_params="k-means++",
                        random_state=0).fit(X_tr)
    ll = g.score_samples(X_te)
    edge = np.linalg.norm(X_te, axis=1, ord=np.inf) > 0.9
    print(f"K={K:3d}  mean log-density {ll.mean():+.4f}"
          f"   interior {ll[~edge].mean():+.4f}   near the edge {ll[edge].mean():+.4f}")

kde = KernelDensity(bandwidth=0.1).fit(X_tr)
ll = kde.score_samples(X_te)
edge = np.linalg.norm(X_te, axis=1, ord=np.inf) > 0.9
print(f"KDE   mean log-density {ll.mean():+.4f}"
      f"   interior {ll[~edge].mean():+.4f}   near the edge {ll[edge].mean():+.4f}")
# true log-density -1.3863
# K=  1  mean log-density -1.7306   interior -1.5358   near the edge -2.5387
# K=  4  mean log-density -1.5749   interior -1.4177   near the edge -2.2268
# K= 16  mean log-density -1.4786   interior -1.3751   near the edge -1.9082
# K= 64  mean log-density -1.4648   interior -1.4192   near the edge -1.6538
# KDE   mean log-density -1.4876   interior -1.4077   near the edge -1.8194

3. Initialization, stopping, and partition stability. Fit the same mixture from fifty starts under two initialization methods. Compare likelihoods and pairwise partition agreement, then continue each method’s best fit with a tighter stopping tolerance.

Rounded likelihood buckets count ranges of scores, not distinct optima. The continuation checks how much the default stopping rule affects the comparison.

import numpy as np
from collections import Counter
from itertools import combinations
from sklearn.mixture import GaussianMixture
from sklearn.metrics import adjusted_rand_score

rng = np.random.default_rng(0)
X = np.vstack([rng.normal([0, 0], 1.0, (400, 2)),
               rng.normal([3, 0], 1.0, (400, 2)),
               rng.normal([1.5, 2.5], 1.0, (400, 2))])
y = np.r_[np.zeros(400), np.ones(400), 2 * np.ones(400)]

for init in ("random", "k-means++"):
    lls, aris, partitions, models = [], [], [], []
    for seed in range(50):
        g = GaussianMixture(3, n_init=1, init_params=init, max_iter=500,
                            random_state=seed).fit(X)
        models.append(g)
        partitions.append(g.predict(X))
        lls.append(g.score(X) * len(X))
        aris.append(adjusted_rand_score(y, g.predict(X)))
    lls, aris = np.array(lls), np.array(aris)
    print(f"{init:10s} loglik  best {lls.max():9.2f}  worst {lls.min():9.2f}"
          f"  spread {lls.max() - lls.min():7.2f}")
    print(f"{'':10s} ARI     best {aris.max():.4f}  worst {aris.min():.4f}"
          f"  median {np.median(aris):.4f}")
    print(f"{'':10s} likelihood buckets (width 0.5 nat) {len(Counter((lls / 0.5).round()))}"
          f"   within 1 nat of best {np.mean(lls > lls.max() - 1):.2f}")
    pairwise = [adjusted_rand_score(a, b) for a, b in combinations(partitions, 2)]
    print(f"           pairwise ARI min {min(pairwise):.4f}, median {np.median(pairwise):.4f}")
    best = models[int(lls.argmax())]
    best.set_params(warm_start=True, tol=1e-8, max_iter=2000)
    best.fit(X)
    print(f"           refined best loglik {best.score(X)*len(X):.2f},"
          f" ARI {adjusted_rand_score(y, best.predict(X)):.4f}, converged {best.converged_}")
# random     loglik  best  -4486.06  worst  -4486.70  spread    0.64
#            ARI     best 0.5986  worst 0.0395  median 0.3566
#            likelihood buckets (width 0.5 nat) 2   within 1 nat of best 1.00
#            pairwise ARI min -0.1409, median 0.2603
#            refined best loglik -4360.97, ARI 0.6398, converged True
# k-means++  loglik  best  -4361.95  worst  -4396.36  spread   34.41
#            ARI     best 0.6761  worst 0.2581  median 0.5714
#            likelihood buckets (width 0.5 nat) 34   within 1 nat of best 0.14
#            pairwise ARI min 0.0249, median 0.5051
#            refined best loglik -4360.97, ARI 0.6398, converged True

References

• Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data via the EM Algorithm. Journal of the Royal Statistical Society, Series B.