Gaussian Mixture Models and Expectation-Maximization

In K-means Clustering every observation belonged to exactly one cluster, and the algorithm committed to that answer the moment it assigned the point. That commitment is convenient, but it hides two things we often care about. The first is doubt: a point sitting between two clusters is not really "80% cluster A and 20% cluster B" to k-means; it is entirely cluster A or entirely cluster B depending on which side of the boundary it happens to fall. The second is shape: k-means measures similarity with a plain Euclidean distance, so its nearest-centroid decision regions have straight boundaries. The squared-distance objective favours compact, roughly isotropic groups, but the regions themselves need not be spherical or equally sized.

A Gaussian mixture model (GMM) replaces that hard machinery with a probability model. We imagine that nature first picked one of $K$ hidden groups for each observation, then drew the observation from a Gaussian distribution belonging to that group. Because the group labels are never observed, they are latent variables. The model does not need to be told the labels; it infers a probability distribution over them, which is exactly the soft assignment that k-means refused to give us.

The price of that richer model is that the likelihood we want to maximize contains the log of a sum, so the general multi-component maximum likelihood equations do not yield a direct closed-form joint parameter estimate like the simpler examples in Maximum Likelihood Estimation. The Expectation-Maximization (EM) algorithm solves this by repeatedly pretending the latent labels are known, taking their expected values under the current parameters, and then re-estimating the parameters as if those expectations were counts. An exact E-step followed by an M-step that increases the auxiliary objective cannot decrease the observed-data likelihood. That guarantee requires the stated updates and valid parameters; the diagonal-addition safeguard used in our practical code does not automatically inherit it. Even exact nondecrease is weaker than finding a global maximum.

By the end of this chapter we will have derived the E-step and M-step from scratch, worked a four-point example entirely by hand, built a numerically stable implementation that handles the underflow problems a naive version would hit, and compared its soft assignments against k-means on clusters that are deliberately overlapping and elliptical.

The most immediate benefit is that the model reports uncertainty. The fitted quantity $\gamma_{ik}$, called the responsibility of component $k$ for observation $i$, is a posterior probability under the fitted mixture: it is the model's chance that this observation was generated by component $k$, conditional on its features and the current parameter estimates. We can threshold it to obtain clusters, keep it as a soft feature for a downstream model, or use it to flag observations that sit between groups and therefore deserve a human look. K-means supplies distances and hard assignments rather than these model-based posterior probabilities. A responsibility is not, by itself, evidence of calibration against real-world groups or a measure of uncertainty in the estimated parameters.

The second benefit is geometric flexibility. Each component carries its own full covariance matrix, so it can be elongated, rotated, and scaled differently from its neighbours. It can also carry its own mixing weight, so one component can be large and diffuse while another is small and tight. Data that k-means would either split in the wrong place or merge into one blob (think of a long diagonal ridge next to a compact ball) is handled naturally. Restricted covariance structures (shared, diagonal, or spherical) let us dial that flexibility back down when data are scarce, which gives us a principled continuum from k-means-like fits to fully general ones.

The third benefit is the fitted density itself. Once fitted, $p(\mathbf{x})$ is a proper probability distribution we can evaluate, sample from, and compare across models. That makes GMMs useful far beyond clustering: for density estimation, for anomaly detection (a point with very low mixture density is unusual), for generative sampling, and as a flexible component inside larger probabilistic models. Because the likelihood is a well-defined criterion, we can select the number of components with information criteria such as BIC rather than eyeballing a scree plot.

For a specified mixture, we can marginalize unobserved coordinates: each component's marginal is Gaussian, and their weighted sum is still a Gaussian mixture. Evaluating this marginal does not require filling in those coordinates. Our code and scikit-learn's GaussianMixture expect complete feature rows; fitting mixtures to incomplete training observations needs an additional missing-data procedure.

The likelihood surface is difficult to optimize. With multiple components the log-likelihood is not concave, so EM can settle near different stationary points. Convergence of parameters requires additional regularity; likelihoods can also approach singular boundary configurations. Different initializations can produce visibly different answers with different likelihoods, and there is no certificate telling us that the best one was found. Multiple restarts and a structured initialization such as k-means are useful safeguards, although neither guarantees the best fit.

For an unrestricted mixture with at least two components, the sample likelihood is unbounded above. Put one component mean on an observation and shrink its covariance to zero while another nondegenerate component gives the remaining observations positive density. The first contribution diverges, even though this is a degenerate density fit. This is not a numerical quirk that better code removes; it is a property of the objective. The standard remedy is to add a small constant $\epsilon$ to the diagonal of every covariance matrix, that is, replace $\boldsymbol{\Sigma}_k$ with $\boldsymbol{\Sigma}_k + \epsilon\mathbf{I}$, where $\epsilon>0$ is a positive regularization strength (often called reg_covar in code) and $\mathbf{I}$ is the $d \times d$ identity matrix. For a positive-semidefinite scatter matrix, this imposes a positive covariance eigenvalue floor and caps the density. It prevents singular factorization but does not establish global optimality or exact EM monotonicity. The size of $\epsilon$ is therefore a real modelling decision, not a rounding detail.

There are further costs. The method assumes each component is Gaussian. Within a fixed feature space, approximating strongly skewed or heavy-tailed groups may require extra Gaussian components, which can make interpretation harder. Components are not the same thing as real-world groups: a single natural group with a non-Gaussian shape may be split across two fitted components, and two natural groups may share one. Component labels are arbitrary: permuting the component indices gives exactly the same likelihood, so comparisons across runs require matching components before comparing parameters. Finally, the computational cost grows with the square of the dimension: evaluating the full-covariance densities costs $O(nKd^2)$ after $O(Kd^3)$ factorizations, and forming the weighted covariance scatters costs $O(nKd^2)$. Here $n$ is the observation count, $d$ the feature count, and $K$ the number of components. For a thousand well-populated clusters in high dimensions, k-means remains far cheaper.

This section derives the mixture density, introduces the latent-variable formulation, and develops the E- and M-step updates.

The derivation starts with the smallest object we need: a single Gaussian bump. In one dimension, with mean $\mu$, positive standard deviation $\sigma>0$ and variance $\sigma^2$, its density is

$$\mathcal{N}(x \mid \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).$$

The fraction in front is not decoration. It is the exact constant that makes the area under the curve equal one, so the result can be treated as a probability density. The exponent measures how far $x$ sits from $\mu$ in units of standard deviations, then squares that distance; the negative sign turns distance into a penalty, so points near the mean get high density and points far away get low density. A single Gaussian therefore encodes one centre and one spread. That is useful, but it cannot describe data with two separated modes.

A mixture removes that limitation by averaging several Gaussian densities together:

$$p(x \mid \boldsymbol{\theta}) = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(x \mid \mu_k, \sigma_k^2),$$

where:

• $\boldsymbol{\theta}$: the full parameter set collecting all means $\mu_k$, variances $\sigma_k^2$, and mixing weights $\pi_k$
• $\pi_k$: the mixing weight of component $k$, satisfying $\pi_k \ge 0$ and $\sum_{k=1}^K \pi_k = 1$, so the weights behave like probabilities
• $\mu_k$: the mean (centre) of component $k$
• $\sigma_k^2$: the variance of component $k$ (one-dimensional case; the multivariate case replaces it with a covariance matrix $\boldsymbol{\Sigma}_k$ below)
• $\mathcal{N}(x \mid \mu_k, \sigma_k^2)$: the Gaussian density of component $k$ evaluated at $x$

The mixing weights describe prior component membership within the model: before seeing an observation's features, $\pi_k$ is its chance of coming from component $k$. These are mixture probabilities, not a Bayesian prior on the fitted parameters. The sum says that the overall density at $x$ is a weighted average of the component densities at $x$: each component contributes in proportion to how much prior mass it owns.

This weighted sum is also a generative recipe. To produce a point, first choose a component $k$ with probability $\pi_k$, then draw the point from that component's Gaussian. The recipe contains a hidden step, and naming that step is what makes inference possible. Let $z_i \in \{1, \dots, K\}$ record which component produced observation $i$. We never observe $z_i$, which is why it is called a latent variable. Writing the model jointly over the observation and its hidden label gives

$$p(\mathbf{x}_i, z_i = k \mid \boldsymbol{\theta}) = \pi_k \, \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k),$$

and the law of total probability recovers the marginal density we use:

$$p(\mathbf{x}_i \mid \boldsymbol{\theta}) = \sum_{k=1}^{K} p(\mathbf{x}_i, z_i = k \mid \boldsymbol{\theta}).$$

That identity is the reason the latent variable earns its place. The marginal likelihood contains a sum over hidden states, and that sum is exactly what blocks a direct closed-form maximum. If we could condition on $z_i$, the sum would disappear; the latent variable lets us separate the difficult marginal problem into a sequence of easier complete-data problems.

Bayes' rule then turns the generative model into an inference procedure. Given the current parameters, the posterior probability that observation $i$ belongs to component $k$ is

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

The denominator is the total density at $\mathbf{x}_i$ contributed by all components, and the numerator is the piece of that total contributed by component $k$. Their ratio is the share of the density that component $k$ explains, which is precisely a probability. These $\gamma_{ik}$ are the soft assignments. They replace the missing labels with graded beliefs instead of forcing each point into one cluster.

Before generalizing, it is worth seeing why the obvious optimization strategy fails. The log-likelihood of $n$ independent observations is

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

We take logs because the likelihood itself is a product of $n$ densities, and products of many small numbers underflow quickly on a computer; logs turn products into sums without moving the maximum, since $\log$ is strictly increasing. The obstacle is the outer logarithm sitting on top of a sum. Differentiating with respect to $\boldsymbol{\mu}_k$ and setting the result to zero produces

$$\sum_{i=1}^{n} \gamma_{ik} \, \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0},$$

because the ratio $\pi_k \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) / \sum_j \pi_j \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)$ is precisely $\gamma_{ik}$. Solving for the mean gives

$$\boldsymbol{\mu}_k = \frac{\sum_{i=1}^{n} \gamma_{ik}\, \mathbf{x}_i}{\sum_{i=1}^{n} \gamma_{ik}}.$$

This looks like a closed-form answer, but it is not: the responsibilities $\gamma_{ik}$ on the right-hand side are themselves functions of the whole parameter set $\boldsymbol{\theta}$, including $\boldsymbol{\mu}_k$. We have a fixed-point equation, not a formula. EM's insight is to freeze the right-hand side at its current value, recompute the left-hand side, and repeat. The E-step supplies the current beliefs about the hidden labels; the M-step turns those beliefs into new parameters.

Now the multivariate case, which is what we will implement. For an observation $\mathbf{x} \in \mathbb{R}^d$ the component density is

$$\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = (2\pi)^{-d/2} \lvert \boldsymbol{\Sigma} \rvert^{-1/2} \exp\!\left(-\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\right),$$

where $\boldsymbol{\mu}$ is a $d$-vector, $\boldsymbol{\Sigma}$ is a $d \times d$ symmetric positive definite covariance matrix, the superscript $\top$ means transpose (turn the column vector into a row vector so the matrix product is defined), $\boldsymbol{\Sigma}^{-1}$ is the matrix inverse, and $\lvert \boldsymbol{\Sigma} \rvert$ is the determinant. The quadratic form $(\mathbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})$ is the squared Mahalanobis distance: it rescales and de-rotates space so that distance is measured in units of the component's own spread. The factor $\lvert \boldsymbol{\Sigma} \rvert^{-1/2}$ does the same job the $1/\sigma$ factor did in one dimension: it divides by the volume the Gaussian occupies, so the density integrates to one.

EM alternates between two problems, each of which is easy on its own. The first problem is inference: given parameters, what do we believe about the hidden labels? The second is estimation: given beliefs about the labels, what parameters are most likely?

Both problems are phrased in terms of the complete-data log-likelihood, the quantity we would maximize if the labels were observed. With $z_{ik}$ a one-hot indicator equal to one when observation $i$ comes from component $k$,

$$\ell_c(\boldsymbol{\theta}) = \sum_{i=1}^{n} \sum_{k=1}^{K} z_{ik} \Big[ \log \pi_k + \log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \Big].$$

where:

• $\ell_c(\boldsymbol{\theta})$: the complete-data log-likelihood, the log-likelihood we would maximize if the latent labels were observed
• $z_{ik}$: a one-hot indicator equal to $1$ when observation $i$ came from component $k$, and $0$ otherwise
• $\pi_k$: the mixing weight of component $k$
• $\mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$: the Gaussian density of component $k$ evaluated at $\mathbf{x}_i$
• $n$: the number of observations, $K$: the number of components

Because the indicators select exactly one term per observation, each component's mean and covariance can be fitted from its assigned observations. The mixing weights still share their sum-to-one constraint and are estimated from the component counts. That is the structure the marginal likelihood destroyed.

Since $z_{ik}$ is unknown, we replace it with its conditional expectation, which is the responsibility: $\mathbb{E}[z_{ik} \mid \mathbf{x}_i, \boldsymbol{\theta}] = \gamma_{ik}$. This produces the auxiliary function

$$Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{t}) = \sum_{i=1}^{n} \sum_{k=1}^{K} \gamma_{ik} \Big[ \log \pi_k + \log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \Big],$$

where $\boldsymbol{\theta}^{t}$ denotes the parameters at iteration $t$. The E-step consists of computing the $\gamma_{ik}$ from the current parameters using Bayes' rule above. Each $\gamma_{ik}$ depends on all components through the denominator. Holding responsibilities fixed separates the mean and covariance estimation by component, while the mixing weights retain their shared normalization.

The M-step maximizes $Q$ over $\boldsymbol{\theta}$. For the means, only the quadratic term depends on $\boldsymbol{\mu}_k$:

$$\log \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) = -\tfrac{d}{2}\log(2\pi) - \tfrac{1}{2}\log\lvert \boldsymbol{\Sigma}_k \rvert - \tfrac{1}{2}(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k).$$

Differentiating the quadratic form with respect to the vector $\boldsymbol{\mu}_k$ uses the chain rule on the squared Mahalanobis distance,

$$\frac{\partial}{\partial \boldsymbol{\mu}_k} (\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = -2\,\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_i - \boldsymbol{\mu}_k),$$

where the minus sign appears because $\boldsymbol{\mu}_k$ enters with a negative sign inside the parenthesis. Setting $\partial Q / \partial \boldsymbol{\mu}_k = \mathbf{0}$ gives

$$\sum_{i=1}^{n} \gamma_{ik} \, \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0}.$$

Because $\boldsymbol{\Sigma}_k^{-1}$ is invertible, we may multiply both sides on the left by $\boldsymbol{\Sigma}_k$ to remove it, leaving $\sum_i \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k) = \mathbf{0}$. Splitting the sum and solving,

$$\boldsymbol{\mu}_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik} \, \mathbf{x}_i, \qquad N_k = \sum_{i=1}^{n} \gamma_{ik}.$$

where:

• $N_k = \sum_{i=1}^{n} \gamma_{ik}$: the effective number of observations assigned to component $k$, formed by summing the responsibilities
• $\gamma_{ik}$: the responsibility of component $k$ for observation $i$, its posterior probability of membership
• $\mathbf{x}_i$: the $i$-th observation

The mean update requires $N_k>0$. When a component has zero effective count, this weighted-mean formula is undefined. We can stop with an error, reinitialize or remove the component, or retain its parameters as an inactive zero-weight component and skip its update. Our implementation below stops with an error. The effective count need not be an integer, because responsibilities are fractional: a point with $\gamma_{ik} = 0.3$ contributes $0.3\mathbf{x}_i$ to the weighted feature sum, then division by $N_k$ converts that sum to the mean. This is the soft-count idea in its simplest form.

For the covariance, we need two matrix-calculus facts. First, $\partial \log\lvert \boldsymbol{\Sigma} \rvert / \partial \boldsymbol{\Sigma} = \boldsymbol{\Sigma}^{-1}$ for a symmetric positive definite matrix, which is the matrix analogue of $\partial \log \sigma^2 / \partial \sigma^2 = 1/\sigma^2$. Second, $\partial (\mathbf{a}^{\top}\boldsymbol{\Sigma}^{-1}\mathbf{a})/\partial \boldsymbol{\Sigma} = -\boldsymbol{\Sigma}^{-1}\mathbf{a}\mathbf{a}^{\top}\boldsymbol{\Sigma}^{-1}$, which follows from differentiating the inverse. Applying both,

$$\frac{\partial Q}{\partial \boldsymbol{\Sigma}_k} = \sum_{i=1}^{n} \gamma_{ik} \left[ -\tfrac{1}{2}\boldsymbol{\Sigma}_k^{-1} + \tfrac{1}{2}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}\boldsymbol{\Sigma}_k^{-1} \right] = \mathbf{0}.$$

Multiply on the left and right by $\boldsymbol{\Sigma}_k$ (which is legal since it is invertible) and use $N_k = \sum_i \gamma_{ik}$:

$$-N_k \boldsymbol{\Sigma}_k + \sum_{i=1}^{n} \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top} = \mathbf{0} \quad\Longrightarrow\quad \boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{i=1}^{n} \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}.$$

where:

• $\boldsymbol{\Sigma}_k$: the $d \times d$ covariance matrix of component $k$
• $\gamma_{ik}$: the responsibility of component $k$ for observation $i$
• $\mathbf{x}_i - \boldsymbol{\mu}_k$: the centered observation, its deviation from component $k$'s mean
• $N_k$: the effective count for component $k$

For this to be a valid nondegenerate Gaussian covariance we need $N_k>0$ and a positive-definite scatter. This is the weighted scatter matrix of component $k$, with each observation's outer product $(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}$ counted $\gamma_{ik}$ times. It is the same estimator we would use with known labels, except that the counts are fractional.

The mixing weights carry a constraint, $\sum_k \pi_k = 1$, so we cannot differentiate and set to zero naively; we add the constraint times a Lagrange multiplier $\lambda$ and differentiate the augmented expression:

$$\mathcal{L} = \sum_{i=1}^{n}\sum_{k=1}^{K} \gamma_{ik} \log \pi_k + \lambda\left(\sum_{k=1}^{K}\pi_k - 1\right), \qquad \frac{\partial \mathcal{L}}{\partial \pi_k} = \frac{N_k}{\pi_k} + \lambda = 0.$$

The first term's derivative is $N_k/\pi_k$ because $\partial \log \pi_k/\partial \pi_k = 1/\pi_k$, and the Lagrange multiplier $\lambda$ is a constant; note that although $\gamma_{ik}$ itself depends on $\pi_k$ through the Bayes formula, in $Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^t)$ the responsibilities are held fixed at their current E-step values, so within this maximization they act as constants on the $\pi_k$ derivative. Rearranging gives $\pi_k = N_k / (-\lambda)$. Summing both sides over $k$ and using $\sum_k \pi_k = 1$ and $\sum_k N_k = \sum_i \sum_k \gamma_{ik} = n$ (each observation's responsibilities sum to one),

$$1 = \sum_{k=1}^{K}\pi_k = \frac{\sum_{k=1}^{K} N_k}{-\lambda} = \frac{n}{-\lambda} \quad\Longrightarrow\quad -\lambda = n \quad\Longrightarrow\quad \pi_k = \frac{N_k}{n}.$$

where:

• $\pi_k$: the mixing weight of component $k$
• $N_k$: the effective count for component $k$
• $n$: the number of observations
• $\lambda$: the Lagrange multiplier enforcing the constraint $\sum_k \pi_k = 1$

So each mixing weight is just the share of the effective count that component $k$ holds.

Collecting everything, one full EM iteration is:

$$\gamma_{ik} = \frac{\pi_k \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_j \pi_j \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)} \qquad\text{(E-step)}$$ $$N_k = \sum_i \gamma_{ik}, \qquad \boldsymbol{\mu}_k = \frac{1}{N_k}\sum_i \gamma_{ik}\mathbf{x}_i, \qquad \boldsymbol{\Sigma}_k = \frac{1}{N_k}\sum_i \gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^{\top}, \qquad \pi_k = \frac{N_k}{n} \qquad\text{(M-step)}$$

The matrix form of the M-step makes the computation obvious. Gather the responsibilities into an $n \times K$ matrix $\mathbf{R}$ whose entry $(i,k)$ is $\gamma_{ik}$, and the data into an $n \times d$ matrix $\mathbf{X}$ whose rows are the observations. Let $\mathbf{1}_n$ be a column vector of ones and $\mathbf{n} = \mathbf{R}^{\top}\mathbf{1}_n$ the column vector of effective counts. Then

$$\mathbf{M} = \operatorname{diag}(\mathbf{n})^{-1} \mathbf{R}^{\top}\mathbf{X}, \qquad \boldsymbol{\Sigma}_k = \frac{1}{n_k} (\mathbf{X} - \mathbf{1}_n \boldsymbol{\mu}_k^{\top})^{\top} \operatorname{diag}(\mathbf{r}_k) (\mathbf{X} - \mathbf{1}_n \boldsymbol{\mu}_k^{\top}),$$

where:

• $\mathbf{R}$: the $n \times K$ responsibility matrix with entry $(i,k) = \gamma_{ik}$
• $\mathbf{X}$: the $n \times d$ data matrix whose rows are the observations
• $\mathbf{1}_n$: a column vector of ones of length $n$
• $\mathbf{n} = \mathbf{R}^{\top}\mathbf{1}_n$: the column vector of effective counts whose $k$-th entry is $N_k$
• $n_k$: the $k$-th entry of $\mathbf{n}$ (same as $N_k$), the effective count of component $k$
• $\mathbf{r}_k$: the $k$-th column of $\mathbf{R}$ (the responsibilities for component $k$)
• $\mathbf{M}$: the $K \times d$ matrix whose $k$-th row is $\boldsymbol{\mu}_k^{\top}$
• $\boldsymbol{\mu}_k$: the mean of component $k$
• $\operatorname{diag}(\cdot)$: turns a vector into a diagonal matrix

Multiplying $\mathbf{R}^{\top}$ by $\mathbf{X}$ sums the rows of $\mathbf{X}$ weighted by the responsibilities (each column accumulates a weighted sum), and dividing by $\mathbf{n}$ turns those weighted sums into weighted averages, giving the $K \times d$ matrix $\mathbf{M}$ whose $k$-th row is $\boldsymbol{\mu}_k^{\top}$. The covariance expression is the same operation applied to the centered data: subtracting $\mathbf{1}_n\boldsymbol{\mu}_k^{\top}$ centers every row, the diagonal matrix reweights the rows, and the product $\mathbf{A}^{\top}\mathbf{W}\mathbf{A}$ accumulates the weighted outer products. In code this is a handful of matrix multiplications rather than $n$ nested loops, which matters because the E-step is repeated until convergence.

The reason EM works at all is that it is coordinate ascent on a lower bound. For a probability distribution $q(z)$ that is positive on the joint model's positive support, multiplying and dividing inside the log by $q(z)$ and applying Jensen's inequality (which says that for a concave function like $\log$, the average of the function is at most the function of the average) gives

$$\log p(\mathbf{x}\mid\boldsymbol{\theta}) = \log \sum_{z} q(z)\frac{p(\mathbf{x}, z \mid \boldsymbol{\theta})}{q(z)} \;\ge\; \sum_{z} q(z) \log \frac{p(\mathbf{x}, z \mid \boldsymbol{\theta})}{q(z)}.$$

The gap between the two sides is the Kullback-Leibler divergence from $q$ to the posterior $p(z \mid \mathbf{x}, \boldsymbol{\theta})$, which is non-negative and equals zero exactly when $q$ is the posterior. The E-step sets $q$ to that posterior, closing the gap so the bound touches the log-likelihood at the current parameters. The M-step then pushes the bound up by maximizing it over $\boldsymbol{\theta}$. Since the log-likelihood is at least the bound, and the bound has increased, the log-likelihood cannot decrease. This proves likelihood nondecrease for the exact bound-increasing updates. It does not prove parameter convergence, a finite maximum, or a global optimum, and it does not apply automatically after an arbitrary covariance modification.

A few properties follow directly from the derivation and shape how the method is used.

• Monotone under exact updates. With valid parameters and an exact E-step followed by a bound-increasing M-step, $\ell(\boldsymbol{\theta}^{t+1}) \ge \ell(\boldsymbol{\theta}^{t})$. The likelihood is non-concave; finite interior limits under suitable regularity are stationary, and they need not be global maxima. Singular boundary behaviour remains possible.
• Unbounded above for unrestricted mixtures with $K\ge2$. A component can collapse onto one observation while another covers the remaining points, driving the sample likelihood to $+\infty$. A single Gaussian whose empirical centered covariance is positive definite does have a finite maximum-likelihood fit; the collapse claim is not a statement about every Gaussian model.
• Label switching. Any permutation of the component indices leaves the likelihood unchanged, so its component labels are arbitrary. Comparing parameter estimates across runs requires matching components first.
• Affine transformation of full-covariance models. For $\mathbf{x}'=A\mathbf{x}+\mathbf{b}$ with invertible $A$, transform each mean to $A\boldsymbol{\mu}_k+\mathbf{b}$ and covariance to $A\boldsymbol{\Sigma}_k A^{\top}$. The density divides by $|\det A|$, so responsibilities remain unchanged and total log-likelihood shifts by $-n\log|\det A|$. This relates corresponding unrestricted fits, not arbitrary numerical runs with fixed regularization. Diagonal models are closed under axis scaling but not general rotations; spherical models permit common uniform scaling and rotations, not arbitrary axis scaling.
• Asymptotic relation to k-means. For fixed observations and centroids, suppose all components have equal fixed weights and common spherical covariance $\sigma^2 \mathbf{I}$. As the externally specified common variance tends to zero, responsibilities at observations with a unique nearest centroid approach nearest-centroid indicators. This varies a prescribed variance across models rather than estimating it in the M-step. For each component with a positive limiting hard-assignment count, its mean update approaches the k-means update. Empty limiting clusters require a separate handling rule. This is the small-variance connection between the mixture and k-means.
• Parameter counts. A full-covariance mixture has $K\big(d + d(d+1)/2\big) + (K-1)$ free parameters: $Kd$ means, $K$ covariance matrices with $d(d+1)/2$ distinct entries each, and $K-1$ free mixing weights. Diagonal covariances cost $2Kd + (K-1)$, spherical ones cost $K(d+1) + (K-1)$, and a tied covariance uses $Kd+d(d+1)/2+(K-1)$ parameters by sharing one covariance across components.

We can carry a full EM iteration by hand on four one-dimensional observations, $x = \{1.0,\ 2.5,\ 4.0,\ 6.0\}$, with two components. To keep the arithmetic transparent we fix the variances at $\sigma_1^2 = \sigma_2^2 = 1.5^2 = 2.25$ and the weights at $\pi_1 = \pi_2 = 0.5$ for the first E-step; the M-step then updates everything.

The normal density at $x$ is $\frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$, where $\mu$ is the mean and $\sigma^2$ is the variance, here $\mu_1 = 2.0$ and $\mu_2 = 5.0$ with $\sigma = 1.5$. With $\sigma = 1.5$ the leading constant is $\frac{1}{1.5\sqrt{2\pi}} \approx 0.26596$. Starting means are $\mu_1 = 2.0$ and $\mu_2 = 5.0$. For $x = 2.5$ the squared distances are $(2.5-2)^2 = 0.25$ and $(2.5-5)^2 = 6.25$; dividing by $2\sigma^2 = 4.5$ gives approximately $0.0556$ and $1.3889$. Evaluating the densities with the unrounded constant and exponents gives approximately $0.25159$ and $0.06632$. The responsibilities follow from the ratio:

$$\gamma_{2,1} = \frac{0.5 \times 0.25159}{0.5 \times 0.25159 + 0.5 \times 0.06632} = \frac{0.25159}{0.31791} = 0.7914, \qquad \gamma_{2,2} = 1 - 0.7914 = 0.2086.$$

Repeating for all four points gives the E-step below. The two outer points are almost certain, while $x = 4.0$ has a mixed posterior assignment, sitting at $\gamma_{3,1} = 0.3392$ (that is, $\gamma_{i1}$ is close to but below one half, so component 2 is slightly favoured). That contested point is useful because a nearest-centroid method still assigns it deterministically, while the mixture retains both possibilities under its fitted probability model.

The log-likelihood at this parameter value is the sum of the logs of the mixture densities. For $x = 1.0$ the mixture density is $\pi_1 \mathcal{N}(1.0 \mid \mu_1, 2.25) + \pi_2 \mathcal{N}(1.0 \mid \mu_2, 2.25) \approx 0.11028133$ using unrounded densities, whose natural logarithm is approximately $-2.20472063$. Summing the four contributions gives $\ell(\boldsymbol{\theta}^{0}) = -8.074$. That is our baseline. The exact M-step cannot lower it; equality is possible at a stationary configuration. In this example we can check the strict increase numerically.

Now the M-step. Effective counts are the column sums: $N_1 = 0.9656 + 0.7914 + 0.3392 + 0.0344 = 2.1306$ and $N_2 = 4 - 2.1306 = 1.8694$. These are not integer counts. They are the total fractional membership each component has earned and provide the denominators for the M-step's weighted averages. The new mean of component 1 is the responsibility-weighted average of the observations,

$$\mu_1 = \frac{0.9656(1.0) + 0.7914(2.5) + 0.3392(4.0) + 0.0344(6.0)}{2.1306} = \frac{4.5073}{2.1306} \approx 2.116,$$

and similarly $\mu_2 = (13.5 - 4.5073)/1.8694 = 4.810$, where we used the fact that the two weighted sums must add to the total $\sum_i x_i = 13.5$. The variances use the same weights applied to squared deviations:

$$\sigma_1^2 = \frac{0.9656(1.0 - 2.116)^2 + 0.7914(2.5-2.116)^2 + 0.3392(4.0-2.116)^2 + 0.0344(6.0-2.116)^2}{2.1306} = 1.428,$$

giving $\sigma_1 = 1.195$, and the same calculation for component 2 gives $\sigma_2^2 = 1.826$, so $\sigma_2 = 1.351$. The weights become $\pi_1 = 2.1306/4 = 0.533$ and $\pi_2 = 0.467$.

The test of the whole procedure is whether the likelihood improved. Recomputing the mixture density at the new parameters gives $\ell(\boldsymbol{\theta}^{1}) = -7.897$, which is larger than $-8.074$. The log-likelihood increases by about $0.177$ in this step, consistent with the nondecrease guarantee for this exact update. The contested point at $x = 4.0$ has also moved: its responsibility for component 2 increases from about $0.661$ to $0.692$. Both means move closer to this point. The posterior change depends on the updated weights and variances as well as the means; distances to the means alone do not determine it.

The cell below reproduces the hand calculation with the same formulas, so we can confirm the arithmetic and inspect the responsibilities directly.

The printed numbers agree with the rounded hand calculation to three decimal places, and the log-likelihood strictly increased. Further exact EM steps cannot lower it, but that alone does not guarantee a finite limiting fit. Our later implementation uses a small-change stopping rule and explicitly checks for a material likelihood decrease. The example therefore demonstrates the full rhythm of EM: infer soft labels, turn them into weighted estimates, and check that the objective moved in the promised direction.

We now build the algorithm from primitives, keeping every step visible. We use NumPy for a Cholesky factorization (to evaluate multivariate normal densities without explicitly forming an inverse) and logsumexp from SciPy (to compute the normalizing constant of the responsibilities in a numerically safe way). The linear algebra functions np.linalg.cholesky and np.linalg.solve are standard NumPy.

Small differences are expected: the two implementations initialize differently, stop at slightly different points, and handle the regularization term in slightly different places. A large gap could reflect different local fits, stopping behaviour, regularization, or an implementation error, which is a reminder to use several restarts whichever library is in play.

The formulas above are mathematically correct and numerically fragile if taken literally. Two things must be changed before they will survive real data.

The first is the responsibility calculation. Forming $\pi_k \mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$ directly and then dividing produces underflow for points far from every component: an exponential of $-800$ is exactly zero in double precision, and dividing zero by zero yields nan. The fix is to compute everything in log space and use the log-sum-exp identity,

$$\begin{aligned} \log \sum_{j=1}^{K} \exp(a_j) &= \log\!\left( e^{m} \sum_{j=1}^{K} \exp(a_j - m) \right) && \text{(collect the common factor)} \\ &= m + \log \sum_{j=1}^{K} \exp(a_j - m) && \text{(use } \log(e^m \cdot s) = m + \log s \text{)} \end{aligned}$$

where $a_j$ are finite real log-joint values and $m = \max_j a_j$.

Factoring out the largest term is a legal algebraic rearrangement: we divide each term by $e^m$ and multiply the total back by $e^m$. It guarantees that the largest exponent is zero, so no term overflows and at least one term is exactly $1$. The responsibility then becomes $\exp(a_k - \text{logsumexp}(\mathbf{a}))$, which is always between zero and one even when the raw densities are unrepresentable. SciPy's logsumexp supplies this stable finite-input calculation. An all-infinite or nan input still needs explicit handling; the identity does not make invalid model parameters meaningful.

The second issue is covariance singularity. A weighted scatter is positive definite when its positive-weight observations span all $d$ feature directions after centring. For a hard-assigned component that requires at least $d+1$ affinely independent points. The responsibility sum $N_k$ is not a rank test: many tiny positive weights can still give a full-rank scatter. Adding $\epsilon\mathbf{I}$ with $\epsilon>0$ bounds the eigenvalues away from zero and makes factorization possible. It also bounds the density of these regularized fits, although adding it after an unconstrained M-step does not guarantee ordinary-likelihood monotonicity.

Initialization deserves as much care as the numerics. Exact bound-increasing updates guarantee likelihood nondecrease, not global optimality. Choose starting values and a restart budget together:

• Random means chosen from the data are simple; some runs land in poor fits.
• K-means partitioning often gives a useful structured start. Fit k-means, use its centroids as the means, its labels to compute initial covariances, and its cluster proportions as initial weights. It can still miss an elongated or overlapping component.
• Multiple restarts with the best final likelihood retained are the pragmatic default; n_init controls the number of starts in scikit-learn.

A useful diagnostic is to track the log-likelihood of every restart: if converged runs have substantially different final values, initialization matters. If some runs have not converged, more iterations may explain part of the spread.

Two further safeguards are worth knowing. Component weights can collapse toward zero when a component becomes responsible for almost nothing, which can reflect a redundant component or a poor initialization but produces $\log \pi_k \to -\infty$ if it reaches exactly zero; flooring the weights at a small positive value keeps the arithmetic finite. And the convergence threshold should be applied to the log-likelihood, not to the parameters, since parameters can drift while the likelihood is flat along a ridge.

The parameters of a Gaussian mixture fit control capacity, shape, and computational effort, and each has a practical consequence.

• n_components ($K$): the number of mixture components. The best attainable fit in nested model families cannot worsen as components are added, but finite EM runs may return worse local fits. Training likelihood alone does not penalize complexity.
• covariance_type: "full" gives each component its own unrestricted matrix ($K d(d+1)/2$ covariance parameters), "tied" shares one matrix across components, "diag" restricts each component to axis-aligned spreads, and "spherical" gives each a single variance. Constraining covariance is the main defence against overfitting when the data are limited, and a shared spherical variance with equal weights has the nearest-centroid hard-assignment rule. Finite-variance soft EM remains a different procedure.
• reg_covar: the $\epsilon$ added to the diagonal. Scikit-learn defaults to $10^{-6}$, the starting value used in this example. Inspect that setting relative to feature scales and degeneracy warnings; a larger value may be needed for an unstable fit.
• init_params and n_init: how the starting point is chosen and how many restarts are attempted. K-means initialization with several restarts is the sensible default.
• max_iter and tol: the iteration cap and the stopping threshold on the log-likelihood. Scikit-learn applies tol to the average lower-bound gain, whereas our teaching implementation checks total log-likelihood change. Inspect the convergence flag and increase the budget if it is exhausted.
• weights_init, means_init, precisions_init: optional user-supplied starting values, useful when there is prior knowledge about where the components should be, or when a specific initialization must be reproduced.

Once fitted, a GaussianMixture object exposes the different ways of using the model, and the distinction between them matters.

• fit(X) estimates the parameters; fit_predict(X) additionally returns the hard labels.
• predict(X) returns the component index with the highest responsibility, which is the hard clustering. It discards the uncertainty that the model went to the trouble of computing.
• predict_proba(X) returns the full responsibility matrix, one row per observation and one column per component. This is usually the more useful output.
• score_samples(X) returns the log mixture density at each observation, which is the natural anomaly score.
• score(X) returns the average log density per observation, the quantity we compared against our own implementation after multiplying by $n$.
• bic(X) and aic(X) return penalized likelihood criteria that can be compared across different values of n_components, since a lower value indicates a better trade-off between fit and complexity.
• sample(n_samples) draws new observations from the fitted density, which is the generative use of the model.

The distinction between predict and predict_proba is central to using the model in practice: the first produces hard GMM labels, which can differ from k-means labels, while the second gives the information that made the mixture model worth fitting.

We begin with the data as generated, coloured by the true component. In a real problem those colours would be unavailable (the whole difficulty is that they are), but seeing them establishes what the model is trying to recover.