Gaussian Mixture Model

A Gaussian mixture model (GMM) describes a probability density as a weighted sum of several Gaussian densities, each with its own mean and covariance. A single Gaussian has one peak and elliptical contours; a mixture can have several peaks and, with enough components, approximate a wide range of shapes. In vision the data are usually pixel colours or local feature vectors, and the model answers two kinds of question: how probable a new colour is under what has been seen (a density model), and which of a few colour populations a pixel probably came from (a soft clustering). Per-pixel mixtures are a classical adaptive background model for video, and colour mixtures for the object and its surroundings drive interactive segmentation such as GrabCut.

This article covers the model, its latent-variable reading, soft assignments, what makes fitting it by maximum likelihood harder than fitting one Gaussian, the covariance types and the choice of the number of components, and its main uses in vision, with a short summary of Expectation–Maximization, the algorithm that fits it.

Definition

A Gaussian mixture model with $K$ components is a probability density on $\mathbb{R}^d$ formed by a convex combination of $K$ Gaussian densities: non-negative mixing weights that sum to one, each multiplying a Gaussian with its own mean vector and covariance matrix. Equivalently, it is the distribution of data produced by a two-stage process: pick one of the $K$ components at random with probabilities given by the weights, then draw the data point from that component’s Gaussian. The identity of the chosen component is not observed; it is a latent (hidden) variable.

Intuition

Take the colours of the pixels in a photograph of a red ball on a lawn under a blue sky. The colours do not form one cloud: they form one cluster near red, one near green and one near blue, each elongated because shading makes each surface brighter or darker along its own colour direction. A single Gaussian fitted to all of them would put its mean at a muddy colour that almost no pixel has and spread its probability over empty regions of colour space. A mixture of three Gaussians can put one ellipsoid on each cluster, with weights proportional to how many pixels each surface covers.

Given such a model and one pixel, Bayes’ rule says how likely each component is to have produced it. A pixel in the middle of the green cluster belongs to the lawn component almost certainly; a pixel whose colour lies where two clusters overlap, say a dark patch of sky that could also be shadowed grass, gets a split verdict such as 60% sky and 40% grass. These posterior probabilities, the responsibilities, are the mixture’s soft assignment; taking the largest gives a hard label, as k-means clustering does.

Formal Definition

The mixture density

The notation follows the site’s convention: bold lowercase column vectors, $p(\cdot)$ for densities, $P(\cdot)$ for probabilities, and a data set $X \in \mathbb{R}^{n \times d}$ with one sample $\mathbf{x}_i^\top$ per row, as scikit-learn stores it. A Gaussian mixture with parameters $\theta = \{\pi_k, \boldsymbol{\mu}_k, \Sigma_k\}_{k=1}^{K}$ has density (Bishop, 2006, §2.3.9)

$$p(\mathbf{x} \mid \theta) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x};\, \boldsymbol{\mu}_k, \Sigma_k), \qquad \pi_k \ge 0, \quad \sum_{k=1}^{K} \pi_k = 1,$$

where each $\Sigma_k$ is symmetric positive definite. The constraints on the mixing weights $\pi_k$ follow from requiring $p$ to be non-negative and to integrate to one, since each component already does. Bishop (2006, §2.3.9) notes that, with enough components and suitable weights, means and covariances, almost any continuous density can be approximated to arbitrary accuracy; Murphy (2022, §3.5.1) states the same for smooth densities.

The latent-variable view

Introduce a discrete random variable $z \in \{1, \dots, K\}$ that records which component generated $\mathbf{x}$ (Murphy, 2022, §3.5; Bishop, 2006, §9.2, uses an equivalent one-of-$K$ binary vector):

$$P(z = k) = \pi_k, \qquad p(\mathbf{x} \mid z = k) = \mathcal{N}(\mathbf{x};\, \boldsymbol{\mu}_k, \Sigma_k).$$

Summing the joint distribution $P(z = k)\,p(\mathbf{x} \mid z = k)$ over $k$ gives back the mixture density, so the two descriptions define the same distribution of $\mathbf{x}$. The latent view gives a sampling procedure, called ancestral sampling: draw $z$ from the categorical distribution with probabilities $\pi_k$, then draw $\mathbf{x}$ from component $z$ (Bishop, 2006, §9.2). Every observed sample $\mathbf{x}_i$ has its own unobserved label $z_i$. Data with the labels attached are called complete, and the observed data without them incomplete (Bishop, 2006, §9.2).

Responsibilities

The posterior probability of each component given an observation follows from Bayes’ rule, with the weights as priors and the component densities as likelihoods (Bishop, 2006, §9.2; Murphy, 2022, §3.5.1):

$$\gamma_{ik} = P(z_i = k \mid \mathbf{x}_i, \theta) = \frac{\pi_k \, \mathcal{N}(\mathbf{x}_i;\, \boldsymbol{\mu}_k, \Sigma_k)}{\sum_{j=1}^{K} \pi_j \, \mathcal{N}(\mathbf{x}_i;\, \boldsymbol{\mu}_j, \Sigma_j)}.$$

Bishop calls $\gamma_{ik}$ the responsibility that component $k$ takes for explaining $\mathbf{x}_i$; Murphy writes $r_{nk}$. Each row of responsibilities sums to one. Using them as fractional memberships is soft clustering; replacing them by the most probable component, $\hat{z}_i = \arg\max_k \gamma_{ik}$, is hard clustering (Murphy, 2022, §3.5.1). The denominator is the mixture density itself, so the same computation that gives the responsibilities also gives the likelihood of each sample.

Properties

Fitting by maximum likelihood

For $n$ independent samples the log-likelihood is (Bishop, 2006, §9.2.1)

$$\ell(\theta) = \sum_{i=1}^{n} \log \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x}_i;\, \boldsymbol{\mu}_k, \Sigma_k).$$

For one Gaussian the logarithm acts directly on the exponential, the log-likelihood is a quadratic in the mean, and the maximum likelihood estimates are the sample mean and the sample covariance with denominator $n$. Here the sum over components sits inside the logarithm, and setting the derivatives to zero gives no closed-form solution (Bishop, 2006, §2.3.9 and §9.2.1). The stationarity conditions do have a suggestive form: each mean must equal a responsibility-weighted average of the data, each covariance a responsibility-weighted covariance about that mean, and each weight the average responsibility of its component. They are not a solution, because the responsibilities depend on the parameters (Bishop, 2006, §9.2.2).

Three further difficulties separate mixtures from a single Gaussian.

• Unbounded likelihood. A component that collapses onto a single data point drives the log-likelihood to infinity, so the global maximum does not exist; Maximum Likelihood Estimation explains why a mixture has this singularity and a single Gaussian does not, and lists the usual remedies. Fitting therefore seeks a well-behaved local maximum. scikit-learn’s reg_covar, below, is one remedy; another is a prior on the covariances, whose MAP estimate Murphy (2022, §8.7.3.4) computes with the same algorithm, and in his Figure 8.26, over five trials with 100 samples each, the MAP fit runs into numerical problems less often than the maximum likelihood fit as the dimension grows.
• Local maxima. Away from the singularities, the log-likelihood generally has several local maxima, and the standard algorithm is not guaranteed to find the largest (Bishop, 2006, §9.2.2). Results depend on the initialization, which is why fitting code runs several initializations or initializes carefully.
• Label switching. Permuting the component labels changes the parameter vector but not the distribution, so every solution comes with $K! - 1$ equivalent ones (Bishop, 2006, §9.2.1). This non-identifiability does not matter for density estimation (Bishop, 2006, §9.2.1), but it does when parameters are compared across fits or across images, or when a posterior over the parameters is sampled, where Murphy (2022, §21.4.1.2) shows the posterior marginals of the means becoming multimodal; an ordering constraint on the means helps in one dimension but has no obvious form in more.

Expectation–Maximization

The standard way to fit a GMM is the Expectation–Maximization (EM) algorithm, for which Bishop (2006, §9.2.2) cites Dempster, Laird and Rubin (1977). It turns the stationarity conditions above into an iteration. The E step computes the responsibilities with the current parameters; the M step re-estimates every component as if the responsibilities were known fractional labels, by the weighted averages just described. Bishop (2006, §9.2.2) states that each E step followed by an M step increases the log-likelihood, and that the iteration is deemed converged when the change in it, or in the parameters, falls below a threshold; EM typically needs many more iterations than k-means, each costlier, so the mixture is commonly initialized from a k-means solution.

Covariance types

The covariance matrices are where most of the parameters are, and constraining them trades flexibility for fewer parameters, faster fitting and stability with little data. scikit-learn’s GaussianMixture offers four types (scikit-learn 1.7.2 documentation), counted here for $K$ components in $d$ dimensions from the library’s own parameter count; every type adds $Kd$ mean parameters and $K - 1$ free weights.

covariance_type	Each component has	Covariance parameters	Contours
"full" (default)	its own general covariance	$K d(d+1)/2$	ellipsoids of any orientation
"tied"	one general covariance shared by all	$d(d+1)/2$	ellipsoids of one shape and orientation
"diag"	its own diagonal covariance	$Kd$	ellipsoids aligned with the axes
"spherical"	its own single variance	$K$	spheres

For RGB pixels ($d = 3$) and $K = 3$, the four types have 29, 17, 20 and 14 free parameters. The choice matters for colour: brightness changes move a surface’s colour along a direction that mixes the channels, so a single surface produces a tilted ellipsoid that a diagonal or spherical component cannot follow. Stauffer and Grimson nevertheless used spherical components, trading accuracy for cost in a model evaluated at every pixel of every frame, and listed full covariances as a future improvement (Stauffer & Grimson, 1999, §2.1 and §5).

The other defaults of GaussianMixture in scikit-learn 1.7.2 are n_components=1, tol=1e-3, reg_covar=1e-6, max_iter=100, n_init=1 and init_params="kmeans". reg_covar is added to the diagonal of every covariance estimate to keep it positive definite; EM stops when the gain in the lower bound, which for this model is the average log-likelihood per sample, falls below tol; and of n_init runs the one with the largest lower bound is kept. With init_params="kmeans" the initial responsibilities are the hard labels of one k-means run, and the alternatives are "k-means++", "random" and "random_from_data" (scikit-learn 1.7.2 documentation and source). When a component’s covariance becomes ill-defined, for instance through singleton or collapsed samples, the library raises an error that suggests fewer components, a larger reg_covar or rescaled data. By our reasoning, the default reg_covar is negligible for 8-bit pixel values, whose noise variances are of order 1 to 100, but not for features scaled to a range of $10^{-3}$, so its effect depends on the units of the data. The method predict_proba returns the responsibilities: its 1.7.2 docstring describes the output as each component’s density, but the code returns the normalized posteriors, which sum to one per sample.

Choosing the number of components

The likelihood on the training data keeps increasing as components are added, so it cannot choose $K$. Penalized-likelihood criteria can. scikit-learn’s bic and aic methods compute, for a fitted model with $p$ free parameters, maximized log-likelihood $\hat{\ell}$ and $n$ samples (scikit-learn 1.7.2 source),

$$\mathrm{BIC} = -2\hat{\ell} + p \ln n, \qquad \mathrm{AIC} = -2\hat{\ell} + 2p,$$

and the model with the lowest value is preferred. These are Murphy’s BIC and AIC losses (2022, §5.2.5): BIC is a simple approximation to the log marginal likelihood, and AIC penalizes complexity less because its penalty does not grow with $n$. Murphy (2022, §21.3.7) contrasts this with k-means: the k-means distortion keeps decreasing as $K$ grows even on validation data, because k-means is a degenerate density of $K$ spikes, while the BIC of a GMM typically falls and then rises. The scikit-learn User Guide (§2.1.1) adds that BIC recovers the true number of components only asymptotically, with much data generated i.i.d. from a Gaussian mixture, names held-out data as the other way to choose it, and points to the variational Bayesian mixture as an alternative that avoids specifying the number.

Relation to k-means

k-means clustering is a limiting case of the mixture. Bishop (2006, §9.3.2) gives every component the covariance $\epsilon I$ with $\epsilon$ fixed and lets $\epsilon \to 0$: the responsibility of the nearest mean tends to one and all others to zero, the EM update of each mean becomes the k-means update, and maximizing the expected complete-data log-likelihood becomes minimizing the k-means distortion. Murphy (2022, §21.4.1.1) states the same result as two approximations, fixed $\Sigma_k = I$ and $\pi_k = 1/K$, and hard instead of soft responsibilities. The differences follow: k-means assigns by Euclidean distance alone, ignoring cluster sizes and shapes, and it estimates no covariances, so it cannot express that one colour population is broad and elongated and another tight. In exchange it is cheaper per iteration and converges in fewer iterations (Bishop, 2006, §9.2.2).

A density model and a clusterer

A fitted mixture is first a density. score_samples gives $\log p(\mathbf{x})$ for any new point, which is what a background model or a colour likelihood in a segmentation energy uses. Reading the components as clusters adds an assumption: that each population of interest is one Gaussian. When a region’s colours are not Gaussian, because of a lighting gradient, texture or pixels that blend two surfaces, the best density may spend several components on one region; when two objects share a colour, one component covers both. A mixture with a good BIC is a good density, not necessarily a segmentation into objects. The second example below shows this.

Examples

Colour populations in a rendered image

The code renders a $96 \times 96$ RGB image with three regions whose colours are drawn from known Gaussians: a blue-grey background, a red disk and a shadow on the background at about 75% of its brightness. Each covariance is elongated along its region’s colour direction (shading) plus isotropic noise. Edges are anti-aliased by supersampling, so pixels on a boundary blend the two regions’ colours in proportion to their area. A three-component mixture is fitted with scikit-learn’s defaults and compared with the truth.

import os
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from sklearn.mixture import GaussianMixture

rng = np.random.default_rng(0)
H = W = 96

# Three regions with known colour Gaussians (RGB, 0-255): a background, a red
# disk and a shadow on the background. Brightness varies along each region's
# colour direction, so each covariance is an elongated ellipsoid.
names = ["background", "red disk", "shadow"]
means = np.array([[90.0, 110.0, 140.0], [190.0, 60.0, 50.0], [67.0, 82.0, 105.0]])
covs = []
for m, s_shade, s_noise in [(means[0], 12.0, 4.0), (means[1], 18.0, 4.0), (means[2], 9.0, 4.0)]:
    a = m / np.linalg.norm(m)
    covs.append(s_shade**2 * np.outer(a, a) + s_noise**2 * np.eye(3))
covs = np.array(covs)

# Coverage of each region per pixel, by 4 x 4 supersampling (anti-aliased edges).
ss = 4
yy, xx = np.mgrid[0:H * ss, 0:W * ss] / ss + 0.5 / ss
disk = (xx - 34.3) ** 2 + (yy - 38.6) ** 2 < 22**2
shadow = (np.abs(xx - 66.2) < 24.6) & (np.abs(yy - 72.7) < 16.1) & ~disk
label_hr = np.where(disk, 1, np.where(shadow, 2, 0))
cover = np.stack([(label_hr == k).reshape(H, ss, W, ss).mean(axis=(1, 3)) for k in range(3)])

# Each pixel blends one colour sample per region by its coverage.
samples = np.stack([rng.multivariate_normal(means[k], covs[k], size=(H, W)) for k in range(3)])
image = np.einsum("khw,khwc->hwc", cover, samples)
X = image.reshape(-1, 3)                       # one pixel per row, as scikit-learn expects
f = cover.reshape(3, -1).T                     # true coverage fractions
pure = f.max(axis=1) == 1.0
truth = f.argmax(axis=1)
print(f"pixels: {len(X)}, mixed boundary pixels: {(~pure).sum()}")

gmm = GaussianMixture(n_components=3, random_state=0).fit(X)   # defaults: full covariances
order = [int(np.argmin(np.linalg.norm(gmm.means_ - m, axis=1))) for m in means]
true_w = f.mean(axis=0)
for k, j in enumerate(order):
    sd_true = np.sqrt(np.linalg.eigvalsh(covs[k]))[::-1]
    sd_fit = np.sqrt(np.linalg.eigvalsh(gmm.covariances_[j]))[::-1]
    print(f"{names[k]:10s} weight {true_w[k]:.3f} -> {gmm.weights_[j]:.3f}  "
          f"mean error {np.linalg.norm(gmm.means_[j] - means[k]):4.1f}  "
          f"axis sds {np.round(sd_true, 1)} -> {np.round(sd_fit, 1)}")
print(f"converged: {gmm.converged_} after {gmm.n_iter_} iterations")

# Soft versus hard assignments on the pure pixels.
R = gmm.predict_proba(X)[:, order]             # responsibilities, columns in true order
hard = R.argmax(axis=1)
conf = R.max(axis=1)
print(f"pure pixels mislabelled by the hard rule: {np.mean(hard[pure] != truth[pure]):.4f}")
for lo, hi in [(0.5, 0.7), (0.7, 0.9), (0.9, 0.99), (0.99, 1.0)]:
    sel = pure & (conf > lo) & (conf <= hi)
    print(f"  max responsibility ({lo:.2f}, {hi:.2f}]: {sel.sum():5d} px, "
          f"mean {conf[sel].mean():.3f}, fraction labelled correctly {np.mean(hard[sel] == truth[sel]):.3f}")
# Mixed pixels on the disk/background edge: responsibility is not area coverage.
edge = ~pure & (f[:, 2] == 0)
r = R[edge, 1]
print(f"disk/background edge pixels: {edge.sum()}, disk coverage between 0.25 and 0.75: "
      f"{np.mean((f[edge, 1] > 0.25) & (f[edge, 1] < 0.75)):.2f}, "
      f"disk responsibility between 0.1 and 0.9: {np.mean((r > 0.1) & (r < 0.9)):.2f}")

# Refit on the pure pixels only, to separate the effect of blended edge pixels.
g2 = GaussianMixture(n_components=3, random_state=0).fit(X[pure])
o2 = [int(np.argmin(np.linalg.norm(g2.means_ - m, axis=1))) for m in means]
print("pure pixels only: weights", np.round(f[pure].mean(axis=0), 3), "->", np.round(g2.weights_[o2], 3),
      " second axis sds", [round(float(np.sqrt(np.linalg.eigvalsh(g2.covariances_[j])[1])), 1) for j in o2])

Output:

pixels: 9216, mixed boundary pixels: 285
background weight 0.664 -> 0.651  mean error  0.2  axis sds [12.6  4.   4. ] -> [12.7  4.1  3.9]
red disk   weight 0.165 -> 0.168  mean error  3.3  axis sds [18.4  4.   4. ] -> [20.2 13.1  4. ]
shadow     weight 0.171 -> 0.181  mean error  2.6  axis sds [9.8 4.  4. ] -> [12.1  4.9  3.8]
converged: True after 8 iterations
pure pixels mislabelled by the hard rule: 0.0102
  max responsibility (0.50, 0.70]:   124 px, mean 0.612, fraction labelled correctly 0.661
  max responsibility (0.70, 0.90]:   239 px, mean 0.826, fraction labelled correctly 0.862
  max responsibility (0.90, 0.99]:   779 px, mean 0.965, fraction labelled correctly 0.981
  max responsibility (0.99, 1.00]:  7789 px, mean 0.999, fraction labelled correctly 1.000
disk/background edge pixels: 124, disk coverage between 0.25 and 0.75: 0.41, disk responsibility between 0.1 and 0.9: 0.02
pure pixels only: weights [0.671 0.163 0.166] -> [0.665 0.163 0.171]  second axis sds [4.0, 4.1, 4.0]

The mixture recovers the three populations from unlabelled pixels: the means are within a few grey levels of the truth, the long axes of the ellipsoids within a few levels too, and the weights within about 0.01. Most of the visible errors come from the 285 blended pixels, which belong to no component of the model. The disk component absorbs blended pixels on its border with the background: its second axis grows from 4 to 13 grey levels to reach towards the background colour; refitted on pure pixels only, every second axis is back to 4.0–4.1 and the weights move closer to the area fractions. A Gaussian mixture of regions is not a model of an image with partial-volume edges, a mismatch that grows with the edge-to-area ratio.

The shadow overlaps the background along the brightness direction, so about 1% of pure pixels fall on the wrong side of the hard decision. The responsibilities say where: among pixels whose largest responsibility is near 0.61 the hard label is right 66% of the time, near 0.83 86% of the time, and near 0.97 98%. When the model is close to the truth the responsibilities behave roughly as calibrated posteriors, and they carry the uncertainty that a hard label throws away. They are not area fractions: 41% of the blended pixels on the disk’s edge are between a quarter and three quarters disk, yet only 2% receive an intermediate responsibility. The likely reason is that a blended colour lies far from both means relative to their narrow spread across the brightness direction, where a small difference in Mahalanobis distance already makes one component far more probable; the disk component’s widened ellipsoid also takes in many of them.

Covariance types and the number of components

Continuing with the same pixels, the code fits 1 to 8 components with each covariance type and compares BIC and AIC, on all pixels and on the pure pixels only, and then repeats the BIC choice for ten random states of the k-means initialization.

# Continues from the previous block (X, pure). BIC and AIC for 1-8 components
# and the four covariance types; lower is better for both.
types = ["full", "tied", "diag", "spherical"]
Ks = range(1, 9)
for data_name, D in [("all pixels", X), ("pure pixels only", X[pure])]:
    print(data_name)
    for ct in types:
        fits = [GaussianMixture(n_components=K, covariance_type=ct, random_state=0).fit(D) for K in Ks]
        bic = np.array([g.bic(D) for g in fits])
        aic = np.array([g.aic(D) for g in fits])
        print(f"  {ct:9s} best K by BIC {Ks[bic.argmin()]}, by AIC {Ks[aic.argmin()]};  "
              f"lowest BIC {bic.min():8.0f};  BIC(K=3) - lowest {bic[2] - bic.min():6.0f}")

# How much does the BIC choice depend on the k-means initialization?
for ct in types:
    picks = []
    for seed in range(10):
        bic = [GaussianMixture(n_components=K, covariance_type=ct, random_state=seed).fit(X).bic(X)
               for K in Ks]
        picks.append(Ks[int(np.argmin(bic))])
    print(f"{ct:9s} BIC choice, all pixels, seeds 0-9: {picks}")

Output:

all pixels
  full      best K by BIC 5, by AIC 5;  lowest BIC   194233;  BIC(K=3) - lowest   3428
  tied      best K by BIC 8, by AIC 8;  lowest BIC   197792;  BIC(K=3) - lowest   5785
  diag      best K by BIC 8, by AIC 8;  lowest BIC   199175;  BIC(K=3) - lowest  11596
  spherical best K by BIC 8, by AIC 8;  lowest BIC   200637;  BIC(K=3) - lowest  12302
pure pixels only
  full      best K by BIC 3, by AIC 3;  lowest BIC   186522;  BIC(K=3) - lowest      0
  tied      best K by BIC 8, by AIC 8;  lowest BIC   187762;  BIC(K=3) - lowest   4893
  diag      best K by BIC 8, by AIC 8;  lowest BIC   189561;  BIC(K=3) - lowest  11263
  spherical best K by BIC 8, by AIC 8;  lowest BIC   190421;  BIC(K=3) - lowest  13246
full      BIC choice, all pixels, seeds 0-9: [5, 5, 5, 5, 5, 5, 5, 5, 4, 5]
tied      BIC choice, all pixels, seeds 0-9: [8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
diag      BIC choice, all pixels, seeds 0-9: [8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
spherical BIC choice, all pixels, seeds 0-9: [8, 8, 8, 8, 8, 8, 8, 8, 8, 8]

With full covariances and only pure pixels, both criteria pick the true three components. Adding the 3% of blended pixels makes BIC prefer five, in nine of the ten initializations, and four in the other: the blended colours are a genuine part of the colour distribution even though they are not objects, and a density that includes them needs more than three full-covariance components. The three constrained types reach their lowest BIC only at the largest $K$ tried, even on pure pixels: they need several axis-aligned or spherical components to tile each tilted ellipsoid, and even then their best BIC is higher than that of three full-covariance components. In this example, choosing the covariance type mattered more than choosing $K$, and the “right” number of components depended on whether the question was the number of surfaces or the shape of the colour density.

Adaptive background subtraction

Background subtraction marks the pixels of a video that differ from a model of the scene without moving objects, a common way to segment moving regions in real time for surveillance and tracking (Stauffer & Grimson, 1999, abstract and §1). A single Gaussian per pixel copes with sensor noise, but a pixel that sees two surfaces in turn, such as swaying branches, water with specular glints or a flickering monitor, has a multimodal history, and a model that does not adapt accumulates errors as lighting and the scene change (§1.1–§2.1). Stauffer and Grimson (1999) modelled the recent history of each pixel as a mixture of $K$ Gaussians, with $K$ from 3 to 5 in their system, each with covariance $\sigma_k^2 I$ to avoid matrix inversions:

$$P(X_t) = \sum_{k=1}^{K} \omega_{k,t}\, \eta(X_t;\, \boldsymbol{\mu}_{k,t}, \Sigma_{k,t}), \qquad \Sigma_{k,t} = \sigma_k^2 I.$$

Running EM over a window of frames at every pixel would be costly, so they used what they call an on-line K-means approximation (§2.1). Each new value $X_t$ is checked against the components until one matches, a match being a value within 2.5 standard deviations. If none matches, the least probable component is replaced by one centred on the new value, with a high initial variance and a low weight. The weights are updated by

$$\omega_{k,t} = (1 - \alpha)\, \omega_{k,t-1} + \alpha M_{k,t},$$

with $M_{k,t} = 1$ for the matched component and 0 for the others, and then renormalized; $1/\alpha$ is the time constant of the adaptation. Only the matched component’s mean and variance change:

$$\boldsymbol{\mu}_t = (1 - \rho)\, \boldsymbol{\mu}_{t-1} + \rho X_t, \qquad \sigma_t^2 = (1 - \rho)\, \sigma_{t-1}^2 + \rho\, (X_t - \boldsymbol{\mu}_t)^\top (X_t - \boldsymbol{\mu}_t), \qquad \rho = \alpha\, \eta(X_t \mid \boldsymbol{\mu}_k, \sigma_k).$$

To decide which components are background, they ranked components by $\omega / \sigma$, which grows as a component gathers evidence and as its variance shrinks, and took the first $B$ of them as the background model (§2.2):

$$B = \arg\min_b \left( \sum_{k=1}^{b} \omega_k > T \right),$$

where $T$ is the minimum fraction of the data the background should account for. A small $T$ usually gives a unimodal background; a larger $T$ lets a repetitive motion, such as leaves or a flag, contribute two or more background colours. A pixel is foreground when it matches none of the background components. Foreground pixels are then grouped by connected components and tracked with a multiple-hypothesis tracker built from Kalman filters. The authors described $\alpha$ and $T$ as the method’s two significant parameters, reported 11 to 13 frames per second for $160 \times 120$ frames on an SGI O2 with an R10000 processor, and reported the system running almost continuously for 16 months on outdoor scenes, in rain and snow (Stauffer & Grimson, 1999, abstract, §1.2 and §3). A side effect of keeping several components is fast recovery: when an object that had become background moves away, the old background component is still in the mixture with a lower weight and is quickly reinstated (§2.1).

The NumPy code implements this update for grey levels on a synthetic video: a static scene with noise of standard deviation 4, a patch that switches at random between two grey levels every frame, and, from frame 250, a bright square that moves slowly across the frame and over the flickering patch. A single running Gaussian per pixel, updated with every frame and with the same $\alpha$ and 2.5-standard-deviation test, is the baseline.

import numpy as np

# Synthetic grey-level video, 64 x 64, 400 frames: a static textured scene with
# sensor noise, a flickering patch that switches at random between two levels
# (like a monitor or leaves against sky), and a bright square crossing the frame
# from frame 250 on, passing over the flickering patch.
rng = np.random.default_rng(1)
H, W, T = 64, 64, 400
yy, xx = np.mgrid[0:H, 0:W]
scene = 60 + 1.5 * xx + 20 * np.sin(yy / 5.0)
flicker = (yy >= 8) & (yy < 24) & (xx >= 36) & (xx < 56)
frames, fg_truth = [], []
for t in range(T):
    f = scene.copy()
    f[flicker] = 80 if rng.random() < 0.5 else 170
    obj = np.zeros((H, W), bool)
    if t >= 250:
        x0 = int(4 + 0.35 * (t - 250))
        obj[10:22, x0:x0 + 12] = True
        f[obj] = 230
    frames.append(f + rng.normal(0, 4, (H, W)))
    fg_truth.append(obj)


def single_gaussian(frames, alpha=0.01, k=2.5, var0=100.0):
    """Running mean and variance per pixel; foreground if beyond k sds."""
    mu, var = frames[0].copy(), np.full(frames[0].shape, var0)
    masks = []
    for x in frames:
        fg = np.abs(x - mu) > k * np.sqrt(var)
        masks.append(fg)
        mu = (1 - alpha) * mu + alpha * x
        var = (1 - alpha) * var + alpha * (x - mu) ** 2
    return masks, np.sqrt(var)


def stauffer_grimson(frames, K=3, alpha=0.01, Tb=0.7, k=2.5, var0=100.0, w0=0.05):
    """Per-pixel mixture with the online update and background rule of
    Stauffer & Grimson (1999), Eqs. 5-9, for grey levels."""
    shape = frames[0].shape
    w = np.zeros((K,) + shape); w[0] = 1.0
    mu = np.zeros((K,) + shape); mu[0] = frames[0]
    var = np.full((K,) + shape, var0)
    masks = []
    for x in frames:
        # Rank components by w / sigma; the first B, whose weights first
        # exceed Tb, form the background model (Eq. 9).
        rank = np.argsort(-w / np.sqrt(var), axis=0)
        w_sorted = np.take_along_axis(w, rank, axis=0)
        B = 1 + (np.cumsum(w_sorted, axis=0) <= Tb).sum(axis=0)
        is_bg = np.zeros((K,) + shape, bool)
        for i in range(K):
            np.put_along_axis(is_bg, rank[i:i + 1], (i < B)[None], axis=0)
        # Match: the first component, in ranked order, within k sds of x.
        within = np.abs(x - mu) <= k * np.sqrt(var)
        within_ranked = np.take_along_axis(within, rank, axis=0)
        any_match = within_ranked.any(axis=0)
        first = np.take_along_axis(rank, within_ranked.argmax(axis=0)[None], axis=0)[0]
        M = np.zeros((K,) + shape, bool)
        np.put_along_axis(M, first[None], any_match[None], axis=0)
        fg = ~(M & is_bg).any(axis=0)
        masks.append(fg)
        # Update weights (Eq. 5) and the matched component (Eqs. 6-8).
        w = (1 - alpha) * w + alpha * M
        eta = np.exp(-0.5 * (x - mu) ** 2 / var) / np.sqrt(2 * np.pi * var)
        rho = alpha * eta
        mu = np.where(M, (1 - rho) * mu + rho * x, mu)
        var = np.where(M, (1 - rho) * var + rho * (x - mu) ** 2, var)
        # No match: replace the least probable component (lowest w / sigma).
        least = rank[-1]
        repl = np.zeros((K,) + shape, bool)
        np.put_along_axis(repl, least[None], (~any_match)[None], axis=0)
        mu = np.where(repl, x, mu); var = np.where(repl, var0, var); w = np.where(repl, w0, w)
        w /= w.sum(axis=0)
    top = w.argmax(axis=0)[None]                 # sd of each pixel's heaviest component
    return masks, np.take_along_axis(np.sqrt(var), top, axis=0)[0]


static = ~flicker
# For each object pixel, how many frames in a row the square has covered it.
dwell, d = [], np.zeros((H, W), int)
for t in range(T):
    d = np.where(fg_truth[t], d + 1, 0)
    dwell.append(d)
for name, (masks, sd) in [("single Gaussian", single_gaussian(frames)),
                    ("mixture, K = 3", stauffer_grimson(frames))]:
    late = range(200, 250)                       # after learning, before the object
    fp_static = np.mean([masks[t][static].mean() for t in late])
    fp_flicker = np.mean([masks[t][flicker].mean() for t in late])
    obj_frames = range(250, T)
    tp = np.mean([masks[t][fg_truth[t]].mean() for t in obj_frames])
    tp_over = np.mean([masks[t][fg_truth[t] & flicker].mean()
                       for t in obj_frames if (fg_truth[t] & flicker).any()])
    short = np.concatenate([masks[t][(dwell[t] > 0) & (dwell[t] < 30)] for t in obj_frames]).mean()
    long_ = np.concatenate([masks[t][dwell[t] >= 30] for t in obj_frames]).mean()
    print(f"{name}: false positives static {fp_static:.4f}, flicker {fp_flicker:.3f}; "
          f"object found {tp:.3f}, over the flicker {tp_over:.3f}\n"
          f"  object pixels covered for < 30 frames found {short:.3f}, for >= 30 frames {long_:.3f}\n"
          f"  median sd at the end: static {np.median(sd[static]):.1f}, flicker {np.median(sd[flicker]):.1f}")

Output:

single Gaussian: false positives static 0.0019, flicker 0.000; object found 0.322, over the flicker 0.001
  object pixels covered for < 30 frames found 0.369, for >= 30 frames 0.000
  median sd at the end: static 4.2, flicker 59.0
mixture, K = 3: false positives static 0.0000, flicker 0.000; object found 0.947, over the flicker 0.907
  object pixels covered for < 30 frames found 1.000, for >= 30 frames 0.579
  median sd at the end: static 9.4, flicker 9.7

Both models are quiet on the static scene and on the flickering patch before the square appears, but for different reasons. The mixture holds the two flicker levels as two background components, whose spread stays near 10. The single Gaussian can only cover both levels with one wide distribution: by our reasoning, two equally frequent levels 90 apart give a standard deviation of about 45, and the square’s own pixels widen it further, to 59 at the end. It is then blind: the square over the flicker is almost never detected. It also misses most of the square elsewhere. By our reasoning, this is because it updates its mean and variance with the square’s pixels too, so frames of a value at least 57 grey levels away inflate its variance until the square falls within 2.5 standard deviations; the mixture instead gives the square a component of its own and leaves the background components unchanged. The mixture’s misses are all pixels that the square has covered for 30 frames or more. While the square covers a pixel, the background component’s weight decays as $(1 - \alpha)^n$, and once it falls below $T$ the square’s component joins the background model; by our reasoning this takes about $\ln T / \ln(1 - \alpha) \approx 35$ frames for $\alpha = 0.01$ and $T = 0.7$, a little less here because the replacement component starts with a weight of 0.05. A slow or stopped object is absorbed into the background after a time set by $\alpha$ and $T$, the behaviour the authors describe for objects that stay long enough to become background. The final standard deviations show a property of the printed update: $\rho = \alpha\eta$ is $\alpha$ times a density value, at most $1/(\sqrt{2\pi}\,\sigma) \approx 0.04$ for $\sigma = 10$, so the matched component’s mean and variance adapt at least 25 times more slowly than the weights. The spread of 10 given to each component at initialization has fallen only to about 9.4 after 400 frames, against a true noise level of 4, so in this example the initial variance acts almost as a fixed threshold. The code checks components in ranked order and replaces the one with the lowest $\omega / \sigma$ when nothing matches, choices the paper leaves open.

Colour models in GrabCut

GrabCut (Rother, Kolmogorov & Blake, 2004) segments an object from a rectangle drawn around it by alternating colour-model estimation with a minimum cut; Graph Cuts describes the energy and the cut. Its colour model is a pair of full-covariance Gaussian mixtures in RGB, one for the foreground and one for the background, typically with $K = 5$ components each, used because colour-space histograms were impractical (§3.1). Each pixel is assigned to a single component of the mixture of its current label, and its data cost is, up to a constant, the negative log of that component’s weighted density, $-\log \pi + \frac{1}{2}\log\det\Sigma$ plus half the squared Mahalanobis distance (Eq. 9). Each iteration assigns components to pixels, re-estimates every component as the sample mean and covariance of its pixels with the weight equal to its share of them, and recomputes the segmentation by a minimum cut; each step minimizes the total energy over one set of variables (component assignments, mixture parameters, segmentation) in turn, so the energy decreases monotonically and the procedure converges at least to a local minimum (§3.2). The assignments are deliberately hard: the authors note that soft assignments would allow EM but cost significantly more computation for negligible practical benefit (§3.1, footnote 1). Their Figure 4 shows the two mixtures overlapping heavily in colour space at initialization and separating as the labelling becomes accurate.

Skin colour and descriptor encodings

Skin detection classifies pixels by their colour likelihood under skin and non-skin models, as in the histogram classifier of Probability and Statistics. Jones and Rehg built both kinds of model from nearly a billion labelled pixels and, comparing histogram and mixture models, found the histograms superior in accuracy and computational cost (Jones & Rehg, 2002, abstract). Mixtures are also fitted to local descriptors rather than pixels: the Fisher vector encoding described in Image Classification represents an image by the gradient of its descriptors’ log-likelihood under a GMM learned from many images.

Common Misconceptions

• “One component per object.” Components model the density, not the scene. In the second example, three regions needed five full-covariance components once edge pixels were included, and eight, the most tried, when covariances were constrained. A region with a lighting gradient or texture needs several components, and two objects of the same colour share one.
• “A responsibility is the fraction of a pixel covered by a component.” It is the posterior probability that the whole pixel came from that component. A blended edge pixel is not described by the model at all, and the first example gave intermediate responsibilities to only 2% of blended pixels.
• “Maximum likelihood gives the best mixture.” The likelihood is unbounded, so its supremum is a degenerate model with a collapsed component; fitting seeks a good local maximum, with regularized covariances or a prior.
• “Soft assignments are always better than hard ones.” They carry uncertainty that hard labels discard, but GrabCut’s authors found hard assignments sufficient at lower cost, and k-means is the hard limit of the same model.
• “Diagonal covariances only lose a little accuracy.” For colour data whose channels move together under shading, the diagonal and spherical types reached their best BIC in the second example only at eight components, the most tried, and even that was worse than three full-covariance components.

Where It Is Used

• Background subtraction: per-pixel adaptive mixtures in the style of Stauffer and Grimson, feeding blob detection and object tracking.
• Segmentation: foreground and background colour models in GrabCut and other graph cut and energy-based methods for image segmentation, usually in an RGB or other colour space.
• Colour models: skin and object colour likelihoods, compared with colour histograms.
• Descriptor encoding: the vocabulary GMM of the Fisher vector in classical image classification.
• Density estimation and clustering: soft clustering of pixels or feature vectors, anomaly scores from $\log p(\mathbf{x})$, with k-means as the usual initialization; fitted by Expectation–Maximization, with k-means as its hard, spherical limit.