Closed-Form Gaussian Mixture M-Step

From Responsibilities to Parameter Updates

In a Gaussian mixture, the component membership of each observation is latent. The EM algorithm handles this by computing posterior responsibilities in its E-step, then maximizing an expected complete-data log-likelihood in its M-step. This tutorial derives the M-step when the component distributions have unrestricted positive-definite covariance matrices. The responsibilities are held fixed throughout this maximization; they act as fractional membership weights, not as labels to be reassigned.

Let $x_1,\ldots,x_n$ be observations in $\mathbb{R}^d$, and suppose the model has $K$ components. Component $k$ has mixing weight $\pi_k$, mean vector $\mu_k$, and covariance matrix $\Sigma_k$. The weights satisfy $\pi_k>0$ and $\sum_{k=1}^K\pi_k=1$; each $\Sigma_k$ is symmetric positive definite. Write $\tau_{ik}$ for the fixed responsibility of component $k$ for observation $i$. Thus $0\leq\tau_{ik}\leq1$ and $\sum_{k=1}^K\tau_{ik}=1$ for each $i$.

Definition: The expected complete-data log-likelihood for fixed responsibilities is $$ Q(\theta\mid\tau)=\sum_{i=1}^n\sum_{k=1}^K\tau_{ik}\left[\log\pi_k+\log\phi_d(x_i;\mu_k,\Sigma_k)\right], $$ where $\theta$ denotes all mixture parameters and $\phi_d(\,\cdot\,;\mu_k,\Sigma_k)$ is the $d$-variate Gaussian density. The M-step maximizes $Q$ over the allowed parameters while keeping the $\tau_{ik}$ fixed.

The relevant weighted sufficient statistics for component $k$ are its effective size $N_k$, weighted sum $s_k$, and weighted second-moment matrix $R_k$. They lead to the weighted mean and centered scatter used in the updates.

N_k=\sum_{i=1}^n\tau_{ik},\qquad s_k=\sum_{i=1}^n\tau_{ik}x_i,\qquad R_k=\sum_{i=1}^n\tau_{ik}x_ix_i^{\mathsf T}.

When $N_k>0$, define $\bar{x}_k=s_k/N_k$ and $S_k=\sum_i\tau_{ik}(x_i-\bar{x}_k)(x_i-\bar{x}_k)^{\mathsf T}$. The identity $S_k=R_k-s_ks_k^{\mathsf T}/N_k$ follows by expanding the centered outer products and using $\sum_i\tau_{ik}x_i=s_k$.

The Closed-Form Updates

Theorem: Suppose every $N_k>0$ and every weighted scatter matrix $S_k$ is positive definite. For unrestricted Gaussian component means and positive-definite covariance matrices, the unique maximizer of $Q$ is $$ \hat{\pi}_k=\frac{N_k}{n},\qquad \hat{\mu}_k=\frac{s_k}{N_k},\qquad \hat{\Sigma}_k=\frac{S_k}{N_k}. $$

The Gaussian log-density is $\log\phi_d(x_i;\mu_k,\Sigma_k)=-\frac d2\log(2\pi)-\frac12\log\det\Sigma_k-\frac12(x_i-\mu_k)^{\mathsf T}\Sigma_k^{-1}(x_i-\mu_k)$. Consequently, the part of $Q$ that depends on component $k$'s weight and Gaussian parameters is

N_k\log\pi_k-\frac{N_k}{2}\log\det\Sigma_k -\frac12\sum_i\tau_{ik}(x_i-\mu_k)^{\mathsf T}\Sigma_k^{-1}(x_i-\mu_k),

up to terms constant in those parameters. The constraints separate: weights are optimized on the probability simplex, while each component's mean and covariance can be optimized using its own weighted observations.

Deriving the Weight and Mean Updates

For the weights, maximize $\sum_kN_k\log\pi_k$ subject to $\sum_k\pi_k=1$ and $\pi_k>0$. Introduce a Lagrange multiplier $\lambda$. The first-order equations and the constraint give

\frac{N_k}{\pi_k}-\lambda=0,\qquad \pi_k=\frac{N_k}{\lambda},\qquad 1=\sum_k\pi_k=\frac{\sum_kN_k}{\lambda}=\frac{n}{\lambda}.

Here $\sum_kN_k=n$, because the responsibilities for each observation sum to one. Therefore $\lambda=n$ and $\hat{\pi}_k=N_k/n$. This is the global maximum: the objective is strictly concave in the positive weights when every $N_k>0$, since its second derivative with respect to $\pi_k$ is $-N_k/\pi_k^2<0$.

For a fixed positive-definite $\Sigma_k$, the mean-dependent term is a weighted quadratic. Differentiating with respect to $\mu_k$ gives

\nabla_{\mu_k}Q =\Sigma_k^{-1}\sum_i\tau_{ik}(x_i-\mu_k) =\Sigma_k^{-1}(s_k-N_k\mu_k).

Setting this gradient to zero yields $\hat{\mu}_k=s_k/N_k$. The quadratic is strictly concave in $\mu_k$, because its Hessian is $-N_k\Sigma_k^{-1}$, which is negative definite. Thus this is the unique maximizing mean for a fixed covariance.

Deriving and Verifying the Covariance Update

After setting $\mu_k=\bar{x}_k$, the covariance-dependent part of $Q$ becomes

q_k(\Sigma_k)=-\frac{N_k}{2}\log\det\Sigma_k -\frac12\operatorname{tr}(\Sigma_k^{-1}S_k).

To verify the maximum globally, set $C_k=S_k/N_k$, which is positive definite by assumption, and write $\Sigma_k=C_k^{1/2}A C_k^{1/2}$ for a positive-definite matrix $A$. Then $\log\det\Sigma_k=\log\det C_k+\log\det A$, and $\operatorname{tr}(\Sigma_k^{-1}S_k)=N_k\operatorname{tr}(A^{-1})$. Comparing the objective at $\Sigma_k$ with its value at $C_k$ gives

q_k(\Sigma_k)-q_k(C_k) =-\frac{N_k}{2}\left[\log\det A+\operatorname{tr}(A^{-1})-d\right].

Let $\lambda_1,\ldots,\lambda_d$ be the positive eigenvalues of $A$. The bracket is $\sum_{j=1}^d(\log\lambda_j+\lambda_j^{-1}-1)$. For $g(\lambda)=\log\lambda+\lambda^{-1}-1$, $g(1)=0$ and $g'(\lambda)=(\lambda-1)/\lambda^2$. Hence $g$ decreases up to $1$, increases after $1$, and is nonnegative, with equality only at $1$. The bracket is therefore nonnegative and vanishes only when $A=I$. This proves that $C_k=S_k/N_k$ is the unique global maximizer of the covariance part, completing the theorem. $\square$

The same update can be calculated from the sufficient statistics without storing each centered outer product:

\hat{\Sigma}_k =\frac{1}{N_k}\left(R_k-\frac{s_ks_k^{\mathsf T}}{N_k}\right).

The divisor is $N_k$, not $N_k-1$. This is the maximum-likelihood update for the weighted Gaussian log-likelihood, not an unbiased sample covariance correction. Responsibilities are fractional weights, and the M-step optimizes the likelihood objective itself.

Worked Updates

Worked Example: A Two-Component Univariate Mixture

Take observations $x=(0,1,4,5)$ and fixed responsibilities for component 1 of $(0.9,0.8,0.2,0.1)$. The component 2 responsibilities are their complements, $(0.1,0.2,0.8,0.9)$. Each observation's two responsibilities sum to one.

For component 1, $N_1=0.9+0.8+0.2+0.1=2$, and $s_1=0.9(0)+0.8(1)+0.2(4)+0.1(5)=2.1$. Thus $\hat{\pi}_1=2/4=0.5$ and $\hat{\mu}_1=2.1/2=1.05$. Its weighted squared deviations sum to $0.9(0-1.05)^2+0.8(1-1.05)^2+0.2(4-1.05)^2+0.1(5-1.05)^2=4.295$. Therefore $\hat{\sigma}_1^2=4.295/2=2.1475$.

For component 2, $N_2=2$ and $s_2=0.1(0)+0.2(1)+0.8(4)+0.9(5)=7.9$, so $\hat{\pi}_2=0.5$ and $\hat{\mu}_2=3.95$. Its weighted squared deviations sum to $0.1(0-3.95)^2+0.2(1-3.95)^2+0.8(4-3.95)^2+0.9(5-3.95)^2=4.295$. Hence $\hat{\sigma}_2^2=2.1475$. Both variances are positive, so these are valid Gaussian covariance updates.

Worked Example: A Full Covariance Update in Two Dimensions

For one component, consider $x_1=(0,0)^{\mathsf T}$, $x_2=(2,0)^{\mathsf T}$, and $x_3=(0,2)^{\mathsf T}$, with responsibilities $1,0.5,0.5$. Then $N=2$, and the weighted sum is $(1,1)^{\mathsf T}$, giving $\hat{\mu}=(0.5,0.5)^{\mathsf T}$.

The centered residuals are $(-0.5,-0.5)^{\mathsf T}$, $(1.5,-0.5)^{\mathsf T}$, and $(-0.5,1.5)^{\mathsf T}$. Weighting the second and third outer products by $0.5$, the scatter is

S= \begin{pmatrix}0.25&0.25\\0.25&0.25\end{pmatrix} +\frac12\begin{pmatrix}2.25&-0.75\\-0.75&0.25\end{pmatrix} +\frac12\begin{pmatrix}0.25&-0.75\\-0.75&2.25\end{pmatrix} = \begin{pmatrix}1.5&-0.5\\-0.5&1.5\end{pmatrix}.

Dividing by $N=2$ gives $\hat{\Sigma}=\begin{pmatrix}0.75&-0.25\\-0.25&0.75\end{pmatrix}$. Its determinant is $0.75^2-0.25^2=0.5>0$, and its first diagonal entry is positive, so this symmetric $2\times2$ matrix is positive definite. The nonzero off-diagonal entries retain the estimated association between the two coordinates.

Worked Example: Fractional Counts and Three Components

For scalar observations $x=(0,1,3,4)$, suppose the fixed responsibility rows are $(0.7,0.2,0.1)$, $(0.2,0.6,0.2)$, $(0.1,0.2,0.7)$, and $(0.4,0.4,0.2)$. Each row sums to one. The effective sizes are $N_1=1.4$, $N_2=1.4$, and $N_3=1.2$, totaling $4$.

The updated mixing weights are $1.4/4=0.35$, $1.4/4=0.35$, and $1.2/4=0.30$. The weighted sums are $s_1=0.2+0.3+1.6=2.1$, $s_2=0.6+0.6+1.6=2.8$, and $s_3=0.2+2.1+0.8=3.1$. Thus the means are $1.5$, $2$, and $31/12$, respectively.

For component 2, the weighted squared deviations from $2$ sum to $0.2(0-2)^2+0.6(1-2)^2+0.2(3-2)^2+0.4(4-2)^2=3.2$, so its variance update is $3.2/1.4=16/7$. This calculation uses soft counts: for example, the first observation contributes $0.2$ to component 2, not a whole observation or none.

Conditions, Interpretation, and a Common Pitfall

The theorem assumes $N_k>0$ and $S_k$ positive definite for every component. A zero effective size leaves that component's mean and covariance unidentified by $Q$; its weight is at the boundary rather than at a positive interior solution. If a weighted scatter is singular, the stated unrestricted positive-definite covariance maximizer does not exist: the objective can increase as covariance eigenvalues in zero-scatter directions approach zero. In finite samples, this can occur when the weighted observations do not span all $d$ dimensions. The closed-form formula still produces a positive-semidefinite scatter divided by $N_k$, but that matrix is not an admissible positive-definite Gaussian covariance in this model.

A common mistake is to round responsibilities to the most likely component before doing the M-step. That replaces the weighted objective by a different, hard-assignment calculation and generally changes every update. The correction is to keep the E-step responsibilities fixed as real-valued weights and use them in $N_k$, $s_k$, and $S_k$. Another frequent error is to divide the scatter by $N_k-1$; the M-step maximizes the Gaussian likelihood, so its covariance update divides by $N_k$.

Key takeaway: With fixed responsibilities, each component behaves like a weighted Gaussian sample: its effective count determines its mixture weight, its weighted mean determines its location, and its weighted centered scatter divided by that count determines its maximum-likelihood covariance.

Check Your Understanding

Use the fixed-responsibility M-step formulas to answer these questions.

Why do the effective component sizes sum to $n$, even when no observation has responsibility exactly zero or one?
For one component, suppose $N_k=5$ and $s_k=(10, -5)^{\mathsf T}$. What is its updated mean?
Why does the unrestricted maximum-likelihood covariance update divide by $N_k$, rather than by $N_k-1$?
What condition on $S_k$ ensures that $S_k/N_k$ is a valid positive-definite covariance update?
What goes wrong if responsibilities are rounded to hard assignments before the M-step?