Tutorials › Multivariate Statistics (Graduate) › Convergence Diagnostics and Restart Strategy

Gaussian Mixtures and EM · Tutorial 642 of 975

Convergence Diagnostics and Restart Strategy

Use complementary diagnostics to distinguish stable, admissible EM fits from runs that merely appear to have converged.

Advanced 12 min read

What You'll Learn

  • Compute relative observed-log-likelihood increments and label-invariant parameter changes
  • Explain why small likelihood increments alone do not establish convergence to a maximum
  • Use covariance eigenvalues to detect nearly degenerate component estimates
  • Compare multiple starts using the same objective and model constraints
  • Match component labels before judging whether candidate fits are equivalent
  • Apply a reproducible rule for retaining and reporting candidate fits

Why EM Needs More Than One Convergence Check

An EM run can make its observed log-likelihood increase at every iteration and still stop at a poor fit. The increments may become small because the algorithm is near a stationary point, because parameters are changing very slowly, or because a component covariance is approaching degeneracy. A single stopping signal cannot distinguish these situations.

For a \(G\)-component Gaussian mixture fitted to observations \(x_1,\ldots,x_n\), write \(\theta=(\pi_1,\ldots,\pi_G,\mu_1,\ldots,\mu_G,\Sigma_1,\ldots,\Sigma_G)\), where the weights sum to one, the means are vectors, and the covariance matrices are positive definite when the model requires nonsingular components. The observed log-likelihood is

$$ \ell(\theta)=\sum_{i=1}^n \log\left\{\sum_{g=1}^G \pi_g\,\phi(x_i;\mu_g,\Sigma_g)\right\}. $$

A practical assessment combines changes in this objective, changes in the fitted parameters, checks on covariance eigenvalues, and evidence from multiple initializations. These checks answer different questions: Did the objective stabilize? Did the fit itself stabilize? Is the estimate numerically admissible? Does another start find a better fit?

Likelihood Increments and Their Limits

Definition: For successive iterates \(\theta^{(t)}\) and \(\theta^{(t+1)}\), the observed log-likelihood increment is \(\Delta_t=\ell(\theta^{(t+1)})-\ell(\theta^{(t)})\). A scale-adjusted increment is \(r_t=|\Delta_t|/(1+|\ell(\theta^{(t)})|)\). The added 1 prevents division by zero, and the denominator makes the diagnostic less sensitive to the overall scale of the log-likelihood.

For exact EM, the increment is nonnegative under the usual conditions below. In computation, a tiny negative increment can arise from rounding; a substantial decrease calls for checking the implementation, the objective being recorded, and whether the update really satisfies the EM or generalized-EM condition. For a generalized EM update, the required condition is that the update does not decrease the relevant \(Q\)-function.

Theorem: Suppose the observed-data density is positive and finite at the observations and at the parameter values under consideration. Let \(Z\) denote the latent allocations, and assume the conditional expectations of the complete-data log-likelihood are finite. If an EM update chooses \(\theta'\) such that \(Q(\theta'\mid\theta)\geq Q(\theta\mid\theta)\), then \(\ell(\theta')\geq\ell(\theta)\).

Proof. For each observation \(x_i\), let \(r_i(z)=P_\theta(Z_i=z\mid X_i=x_i)\) be the current posterior allocation probabilities. For any candidate parameter \(\vartheta\), the identity for a discrete latent allocation is

$$ \log p_{\vartheta}(x_i) = \sum_z r_i(z)\log\frac{p_{\vartheta}(x_i,z)}{r_i(z)} + \sum_z r_i(z)\log\frac{r_i(z)}{P_{\vartheta}(Z_i=z\mid X_i=x_i)}. $$

The second sum is the Kullback–Leibler divergence from \(r_i\) to the candidate posterior, so it is nonnegative. Sum over \(i\) and subtract the identity at \(\vartheta=\theta\). At \(\vartheta=\theta\), the candidate posterior equals \(r_i\), making every divergence zero. The result is

$$ \ell(\theta')-\ell(\theta) = Q(\theta'\mid\theta)-Q(\theta\mid\theta) + \sum_{i=1}^n \operatorname{KL}\!\left( P_\theta(Z_i\mid x_i)\,\middle\|\,P_{\theta'}(Z_i\mid x_i) \right) \geq 0. $$

The first difference is nonnegative by the update condition, and each divergence is nonnegative. This proves the claim. \(\square\)

Monotonicity is a useful check, not a convergence certificate. Nor is there a universal threshold for \(r_t\): its practical value depends on numerical precision, data scale, and the required accuracy. A defensible stopping rule specifies a tolerance in advance and asks for stability over several consecutive iterations rather than accepting one small increment.

Worked Example: Checking Relative Likelihood Increments

Suppose a run records \(\ell(\theta^{(t)})=-1200.00\), then \(-1199.20\), then \(-1199.19\), then \(-1199.1895\). Use \(r_t=|\Delta_t|/(1+|\ell(\theta^{(t)})|)\) and a stated threshold of \(10^{-5}\).

The three increments are \(0.80\), \(0.01\), and \(0.0005\), respectively. For the second increment, \(r_t=0.01/(1+1199.20)=0.01/1200.20\approx 8.33\times10^{-6}\). For the third, \(r_t=0.0005/(1+1199.19)=0.0005/1200.19\approx4.17\times10^{-7}\). Both are below \(10^{-5}\), while the first increment is not. Thus the likelihood-increment criterion is met for the final two steps, but this alone does not justify declaring the fit converged; parameter stability and covariance checks remain necessary.

Parameter Stability and Label Matching

A parameter-change diagnostic asks whether the fitted mixture has stopped moving, not just whether its objective has changed little. Raw Euclidean distances can be misleading because weights, mean coordinates, and covariance entries have different scales. One option is a scaled, label-invariant distance. Choose positive scales \(s_\pi,s_\mu,s_\Sigma\) in advance, and let \(\mathfrak S_G\) be the set of permutations of \(G\) component labels. Define

$$ d_s^2(\theta,\theta') = \min_{\tau\in\mathfrak S_G} \sum_{g=1}^G \left[ \frac{(\pi_g-\pi'_{\tau(g)})^2}{s_\pi^2} + \frac{\|\mu_g-\mu'_{\tau(g)}\|_2^2}{s_\mu^2} + \frac{\|\Sigma_g-\Sigma'_{\tau(g)}\|_F^2}{s_\Sigma^2} \right]. $$

Here \(\|\cdot\|_2\) is the Euclidean norm and \(\|\cdot\|_F\) is the Frobenius norm. The minimization pairs components in the most comparable way. Scales should reflect meaningful changes for the application; alternatively, variables can be standardized before fitting and the resulting conventions stated. The distance is a diagnostic, not a universal measure of practical equivalence.

Result: If \(\theta'\) is exactly the same mixture as \(\theta\) with only its component labels permuted, then \(d_s(\theta,\theta')=0\). Conversely, if \(d_s(\theta,\theta')=0\), the parameter tuples agree exactly after some component permutation.

Justification. An exact relabeling supplies a permutation for which every paired weight, mean, and covariance agrees, so the sum is zero and its minimum is zero. In the reverse direction, the minimum is over a finite set of permutations. A zero minimum therefore occurs for at least one permutation; since all scales are positive, every nonnegative squared term in that sum must be zero. Each paired parameter is equal, which is precisely equality up to relabeling. \(\square\)

Worked Example: Recognizing an Exact Relabeling

Consider two two-component fits in one dimension. Fit A has weights \((0.4,0.6)\), means \((0,10)\), and variances \((1,4)\). Fit B, in its displayed order, has weights \((0.6,0.4)\), means \((10,0)\), and variances \((4,1)\).

Match A's first component, with mean 0, to B's second component, also with mean 0. Their weights are both \(0.4\), and their variances are both 1. Match A's second component, with mean 10, to B's first component, also with mean 10. Their weights are both \(0.6\), and their variances are both 4. This permutation makes every term in the distance zero. The fits are identical up to relabeling; comparing components only in displayed order would incorrectly suggest large parameter changes.

Covariance Eigenvalues as an Admissibility Check

For a symmetric positive-definite covariance estimate \(\Sigma_g\), all eigenvalues are positive. Its smallest eigenvalue measures variance in the least-variable direction; an extremely small value can signal a nearly singular component covariance and unstable density calculations. In a mixture, a component can become concentrated around a small subset of observations, producing a very high local density. This is why likelihood improvement and covariance health must be checked together.

An eigenvalue threshold must be interpreted in the measurement scale. If variables have very different units, a numerical floor applied to an unscaled covariance matrix can be arbitrary. Standardize the variables or set a scale-aware lower bound, and be explicit about the constraint. When estimation imposes \(\lambda_{\min}(\Sigma_g)\geq\epsilon\), a candidate violating that constraint is inadmissible, even if its computed likelihood is attractive. A candidate at the floor is feasible but should be flagged: the constraint may be materially shaping the estimate.

Worked Example: Finding a Nearly Singular Covariance

Suppose a fitted two-dimensional component has covariance estimate \(\widehat\Sigma=\begin{pmatrix}0.09&0.06\\0.06&0.04\end{pmatrix}\). Its trace is \(0.13\) and its determinant is \(0.09(0.04)-0.06^2=0\). For a symmetric \(2\times2\) matrix, the eigenvalues have sum equal to the trace and product equal to the determinant. They are therefore \(0.13\) and \(0\).

This matrix is singular, not positive definite, so it is not a valid nonsingular Gaussian covariance. If the intended constraint is \(\lambda_{\min}\geq10^{-4}\), the estimate violates it because \(0<10^{-4}\). In practice, floating-point calculations might return a tiny positive eigenvalue instead of exact zero; it would still fail the stated floor. The diagnostic should use eigenvalues computed from the symmetric part of the numerical estimate and should record the minimum value, not merely whether the software returned a density.

Multiple Starts and a Retention Rule

EM can converge to different stationary fits from different initial values. Multiple starts, such as those based on moment or clustering methods, probe this sensitivity; no finite collection of starts proves that the global maximum has been found. The comparison is meaningful only when every run uses the same data, component count, constraints, and objective. If a penalty is part of the estimation target, compare the same penalized objective across runs and report the observed log-likelihood separately.

A reproducible retention rule separates admissibility from ranking. First reject runs that fail the specified stopping checks, produce nonfinite objective values, or violate the covariance constraints. Next rank the remaining candidates by the common objective being optimized. Retain the highest-ranked admissible fit, while keeping other admissible candidates for sensitivity checks. If several distinct fits have very similar objective values, report that sensitivity rather than treating a small numerical difference as decisive.

Worked Example: Selecting Among Starts

Suppose four runs use the same two-component model with covariance floor \(10^{-4}\). Runs A and B meet the likelihood and parameter-stability criteria, have minimum covariance eigenvalues above the floor, and both have log-likelihood \(-100.4\). Their parameter estimates are exact label permutations as in the preceding example. Run C has log-likelihood \(-99.8\) but a minimum eigenvalue of \(10^{-6}\), so it violates the floor. Run D is admissible and converged, with log-likelihood \(-100.7\).

A and B represent the same fit after label matching, and their objective values agree. C has the largest numerical log-likelihood, but it is not an admissible candidate because \(10^{-6}<10^{-4}\). D is admissible but ranks below A and B because \(-100.7<-100.4\). The retained fit is therefore A (or its equivalent relabeling B), with C recorded as a rejected, constraint-violating run. This conclusion follows from the stated rule, not from likelihood ranking alone.

Common Pitfall: Treating a Small Increment as Proof of a Maximum

A tiny or zero likelihood increment does not establish either parameter convergence or maximality. For example, let \(f(\theta)=-(\theta-1)^2\). The constant sequence \(\theta_t=0\) has objective \(f(0)=-1\) at every step, so every increment is zero; nevertheless, the maximum is \(f(1)=0\), attained at \(\theta=1\). Thus zero increments do not imply that a maximum has been reached.

Nor do increments alone imply parameter convergence. For the same objective, the sequence alternating between \(\theta_t=0\) and \(\theta_t=2\) has \(f(\theta_t)=-1\) at every iteration, so all objective increments are zero, while the parameter sequence does not converge. These examples are logical cautions about the diagnostic: the objective can be uninformative about movement or optimality. In an actual mixture analysis, use objective increments alongside label-matched parameter changes, eigenvalue checks, and multiple starts. A stable run is evidence of a reliable numerical fit under the chosen constraints, not a proof of a global optimum.

Key takeaway: Retain a candidate only after checking stability in both objective and parameters, verifying covariance admissibility, and comparing it with other starts under the same objective. Match labels before comparing component parameters, and state thresholds and constraints so the decision can be reproduced.

Check Your Understanding

Answer these questions using the diagnostic principles developed above.

  1. Why is a scale-adjusted likelihood increment useful, and why does a small value not by itself certify a maximum?
  2. What assumptions make the EM likelihood monotonicity argument valid, and what condition must a generalized EM update satisfy?
  3. How can two fits have zero label-matched parameter distance while their displayed component parameters differ?
  4. A covariance estimate has minimum eigenvalue \(0.00002\), while the specified floor is \(0.0001\). Should its likelihood determine the retained fit? Explain.
  5. What should be held fixed when ranking fits obtained from multiple starts?
Clinical Biostats