Gaussian Mixture Models and Expectation-Maximization › Practical Implications

Gaussian mixtures are a reasonable choice when observations are continuous and the groups you expect are roughly elliptical (elongated, rotated, or differently scaled), and when graded membership is more useful than a single label. For a segmentation task we can retain overlapping membership probabilities. The fitted density also supplies log-density scores for anomaly detection, supports generative sampling, and can serve as a density component in a larger probabilistic model. Full covariance is the version that can use correlations between features and can let each component have its own orientation and aspect ratio, which is the flexibility k-means lacks. When the scientific question only needs a partition, Euclidean feature distance is meaningful, and the clusters are compact, roughly spherical, and of comparable spread, k-means remains the better default.

Cost and sample size determine whether that flexibility is affordable. A GMM iteration costs $O(nKd^2)$ plus $O(Kd^3)$ for covariance factorizations, compared with $O(nKd)$ for k-means, and a full-covariance model has $Kd$ means and $Kd(d+1)/2$ covariance entries to estimate. Each component therefore needs a meaningful number of effective observations: a small $N_k$ is a warning about statistical stability, but nonsingularity depends on the affine span of positive-weight observations, not on $N_k>d+1$. Full covariance estimation can still be unstable with few well-supported observations, especially under overlap or noisy features. When counts fall short, constrain the covariance, reduce dimensionality before fitting (the subject of Dimensionality Reduction), or reduce $K$. Large $n$ and high $d$ with nearly spherical groups are the regime where k-means is usually preferable. The cost of a full-covariance fit also grows quickly with dimension, so a full model in dozens of dimensions should be justified by both the sample size and evidence that correlated, non-spherical components are present.

The data requirements concern the conditional distributions within components. Features should be continuous or transformed to a scale on which a Gaussian is a plausible local approximation; strong skew, heavy tails, or discrete counts are better handled by transforming the features, adding components, or choosing a different model. Outliers can pull means and inflate covariances, and near-duplicates can encourage degenerate components, so inspect fitted responsibilities and effective counts rather than trusting a single likelihood value. Scaling has a precise role. For corresponding full-covariance parameters under an invertible affine transformation, responsibilities are unchanged but density values change by the Jacobian factor. Diagonal models remain closed under feature-wise rescaling; spherical models generally do not. Fixed isotropic regularization, initialization and stopping can change even a full-covariance numerical fit. Scaling is therefore a modelling and numerical choice. Fit any learned preprocessing on training observations and apply it unchanged to held-out observations; compare density scores only in the same coordinates.

A practical workflow starts from a structured initialization and several restarts. K-means centroids are one useful starting point, and trying n_init between 5 and 10 is a practical initial budget rather than a guarantee; keep the run with the highest final log-likelihood. Our example uses reg_covar=1e-6, allows 300 iterations, and uses a total-log-likelihood tolerance of 1e-6. Other data and the scikit-learn average-gain convention may require different settings. After fitting, check that:

• no component has a vanishing weight or an effective count near zero;
• no covariance eigenvalue is sitting on the regularization floor;
• the log-likelihood history is monotone and has plateaued;
• the responsibilities reflect the expected component overlap; nearly hard assignments can be appropriate for well-separated components;
• independent restarts produce similar parameters after matching components.

These checks are diagnostic, not proof of correctness: EM can converge cleanly to a poor local optimum, and a well-behaved finite history does not establish stationarity or global optimality.

Deployment adds concerns that are easy to overlook during exploratory fitting. Store the fitted weights, means, covariances, the preprocessing used for the features, and the random seed, because a later refit without them may produce a different labelling of the components. Monitor held-out average log density rather than training likelihood, and refit or warm-start when the distribution drifts. Component labels are not stable across fits (permuting them leaves the likelihood unchanged), so compare components only after matching them by location, shape, or predictive role. If the model is used for anomaly detection, set the threshold from a validation sample and recheck it when the model is updated. Missing values deserve a note of caution: the Gaussian marginal over observed coordinates is well defined, but not every software implementation accepts incomplete rows, so a separate imputation or marginalization step may be required.

The covariance type is the first structural decision because it controls both flexibility and the number of parameters that must be estimated. Full covariance is the most flexible and the most expensive; it can represent any elliptical shape, at $d(d+1)/2$ parameters per component. Tied covariance shares a single shape across all components, which imposes the same covariance size, shape and orientation on every component while allowing different means and weights, and it reduces the number of covariance parameters. Diagonal covariance restricts ellipses to be axis-aligned, a strong assumption that a global principal component rotation does not automatically justify: removing marginal correlation need not remove correlation within every component. Spherical covariance gives every component a single variance and is the closest mixture analogue to k-means. Choose the most constrained structure that captures the visible structure, and relax it only when the constrained fit leaves systematic correlation or curvature in the residuals or when held-out density clearly improves.

Choosing the number of components is a model selection problem, and the training likelihood cannot do it: the best attainable likelihood within nested model families cannot decrease when components are added, but a finite EM fit can return a worse local solution. The standard criteria penalize complexity. The Bayesian information criterion is

$$\mathrm{BIC} = -2\ell(\hat{\boldsymbol{\theta}}) + p \log n,$$

where:

• $\ell(\hat{\boldsymbol{\theta}})$: the maximized log-likelihood when a finite maximum-likelihood fit exists
• $p$: the number of free parameters
• $n$: the number of observations

The numerical comparison below evaluates this expression at the finite parameters returned by EM; it does not certify that the likelihood was maximized. The penalty per parameter grows with sample size. The Akaike criterion replaces $\log n$ with $2$. Here $n=400$, so BIC penalizes each parameter more heavily because $\log 400>2$; that ordering does not hold for every sample size. These conventional criteria are useful model-selection heuristics, and their regular asymptotic justification needs care for singular mixture models; neither is a guarantee. A practical approach is to use these criteria as a guide across a manageable grid of plausible candidate values. Cross-validated log-likelihood is a reasonable alternative, provided it is computed on held-out data with the density evaluated there rather than on the training set. Because each candidate fit can land in a different local optimum, compare models only after giving each $K$ several restarts and retaining its best likelihood.