Tutorials › Multivariate Statistics (Graduate) › Covariance Constraints and Penalized Mixture Estimation

Gaussian Mixtures and EM · Tutorial 640 of 975

Covariance Constraints and Penalized Mixture Estimation

Learn how covariance constraints and penalties prevent Gaussian mixture components from collapsing, and how each approach changes the fitted covariances.

Advanced 11 min read

What You'll Learn

  • Formulate an eigenvalue-floor constraint for each Gaussian mixture covariance.
  • Derive the constrained covariance update by truncating empirical eigenvalues.
  • Distinguish a lower bound on a structured component’s scalar scale from a general eigenvalue floor.
  • Derive a covariance update under an inverse-trace penalty.
  • Compare how floors and penalties affect the fitted covariance and the meaning of the estimate.
  • Recognize what these regularizations do not guarantee about mixture fitting.

Why Constrain Covariances?

In a finite Gaussian mixture, a component can concentrate around one or a few observations by shrinking one or more covariance eigenvalues. Its density at those observations can then become arbitrarily large. As established in the discussion of unbounded Gaussian mixture likelihoods, the ordinary likelihood may therefore have no finite maximum. Covariance constraints and penalties address this instability, but they do so in different ways: a hard constraint excludes part of the parameter space, while a penalty changes the objective being optimized.

Let \(x_1,\ldots,x_n\in\mathbb{R}^d\) be observations and consider a mixture with \(K\) components, weights \(\pi_k\), means \(\mu_k\), and positive-definite covariance matrices \(\Sigma_k\). The weights satisfy \(\pi_k\geq 0\) and \(\sum_{k=1}^K\pi_k=1\). In an EM update, let \(r_{ik}\) denote the current responsibility of component \(k\) for observation \(i\). Define its effective membership and weighted scatter by

$$ N_k=\sum_{i=1}^n r_{ik}, \qquad C_k=\sum_{i=1}^n r_{ik}(x_i-\mu_k)(x_i-\mu_k)^\mathsf{T}. $$

For fixed responsibilities and means, the covariance-dependent part of the expected complete-data log-likelihood is

$$ Q_k(\Sigma_k) =-\frac{N_k}{2}\log|\Sigma_k| -\frac{1}{2}\operatorname{tr}(\Sigma_k^{-1}C_k) +\text{constant}. $$

When \(N_k>0\), the unconstrained update is \(C_k/N_k\). If that scatter is singular, or nearly so, the update may be singular or numerically unstable. The choices below alter this update in distinct, interpretable ways.

An Eigenvalue Floor

Definition: An eigenvalue floor with threshold \(\varepsilon>0\) restricts each component covariance to satisfy \(\Sigma_k\succeq\varepsilon I_d\), where \(I_d\) is the \(d\times d\) identity matrix. Equivalently, every eigenvalue of \(\Sigma_k\) must be at least \(\varepsilon\).

This is a hard constraint on every direction in the data space. It does not require a spherical covariance: eigenvalues may differ and eigenvectors may point in any direction, provided no eigenvalue falls below the threshold.

Theorem: Suppose \(N_k>0\), \(C_k\) is positive semidefinite, and \(\varepsilon>0\). Let \(C_k/N_k=U\operatorname{diag}(s_1,\ldots,s_d)U^\mathsf{T}\), with \(U\) orthogonal and \(s_j\geq0\). A maximizer of \(Q_k(\Sigma_k)\) subject to \(\Sigma_k\succeq\varepsilon I_d\) is $$ \widehat{\Sigma}_{k,\varepsilon} = U\operatorname{diag}\bigl(\max(s_1,\varepsilon),\ldots,\max(s_d,\varepsilon)\bigr)U^\mathsf{T}. $$

Proof. Write \(S=C_k/N_k\). Apart from a positive multiplicative constant and an additive constant, maximizing \(Q_k\) is equivalent to minimizing \(\log|\Sigma|+\operatorname{tr}(\Sigma^{-1}S)\) over \(\Sigma\succeq\varepsilon I_d\). For fixed eigenvalues of \(\Sigma\), the trace term is minimized by aligning its eigenvectors with those of \(S\), pairing the largest eigenvalue of \(S\) with the largest eigenvalue of \(\Sigma\). This follows from the trace rearrangement inequality applied to \(S\) and \(\Sigma^{-1}\): the inverse eigenvalues are in the opposite order. Thus it suffices to minimize separately, for each \(s_j\), the scalar function \[ g_j(\lambda)=\log\lambda+\frac{s_j}{\lambda}, \qquad \lambda\geq\varepsilon. \] Its derivative is \(g_j'(\lambda)=(\lambda-s_j)/\lambda^2\). If \(s_j\geq\varepsilon\), the function decreases up to \(s_j\) and increases afterward, so its minimum is at \(s_j\). If \(s_j<\varepsilon\), it is increasing throughout the allowed range, so its minimum is at \(\varepsilon\). This gives \(\lambda_j=\max(s_j,\varepsilon)\) in every direction and proves the result. \(\square\)

Worked Example: Flooring a Singular Empirical Covariance

Suppose a component has \(N_k=5\) and weighted scatter \(C_k=\operatorname{diag}(20,0.2)\). Then \(C_k/N_k=\operatorname{diag}(4,0.04)\). With an eigenvalue floor \(\varepsilon=0.25\), the two empirical eigenvalues are updated as

$$ \max(4,0.25)=4, \qquad \max(0.04,0.25)=0.25. $$

The constrained estimate is therefore \(\operatorname{diag}(4,0.25)\). The direction with substantial observed spread is unchanged; the direction with very little spread is raised to the specified minimum. The resulting covariance is positive definite even though the empirical covariance was close to singular.

A Lower Bound on a Structured Component Scale

A scale bound is not automatically the same as an unrestricted eigenvalue floor. To see the distinction, suppose the covariance has a fixed shape matrix \(A_k\), which is positive definite, and a scalar scale parameter \(\tau_k^2\): \(\Sigma_k=\tau_k^2 A_k\). Impose \(\tau_k^2\geq a\), where \(a>0\). This restricts one scalar parameter within the structured model; it does not let the data freely adjust each covariance eigenvalue.

Result: For \(N_k>0\), the maximizer of \(Q_k\) under \(\Sigma_k=\tau_k^2 A_k\) and \(\tau_k^2\geq a\) is $$ \widehat{\tau}_k^2 = \max\left( a,\, \frac{\operatorname{tr}(A_k^{-1}C_k)}{dN_k} \right). $$

Derivation. Substituting \(\Sigma_k=\tau_k^2 A_k\) into \(Q_k\), the terms depending on \(\tau_k^2\) are \[ -\frac{N_kd}{2}\log(\tau_k^2) -\frac{1}{2\tau_k^2}\operatorname{tr}(A_k^{-1}C_k). \] Set \(b=\operatorname{tr}(A_k^{-1}C_k)/(dN_k)\). Up to a positive factor and a constant, maximizing this expression is equivalent to minimizing \(h(t)=\log t+b/t\) over \(t\geq a\). Since \(h'(t)=(t-b)/t^2\), its minimum is at \(b\) if \(b\geq a\), and at the boundary \(a\) otherwise. Hence the stated update follows. \(\square\)

Worked Example: A Bound on Scale with Fixed Shape

Take dimension \(d=2\), effective membership \(N_k=4\), and fixed shape \(A_k=\operatorname{diag}(2,1/2)\). Let \(C_k=\operatorname{diag}(2,2)\) and require \(\tau_k^2\geq0.8\). First calculate \[ A_k^{-1}=\operatorname{diag}(1/2,2), \qquad \operatorname{tr}(A_k^{-1}C_k)=1+4=5. \]

The unconstrained scale estimate is \(5/(2\cdot4)=0.625\), which is below the lower bound. Thus \(\widehat{\tau}_k^2=0.8\), and the fitted covariance is \(0.8A_k=\operatorname{diag}(1.6,0.4)\). This update preserves the prescribed shape ratio: the larger variance is always four times the smaller one. An unrestricted eigenvalue-floor update would not impose that ratio.

In a fixed-shape model, the scale bound constrains spread through the chosen parameterization. If the shape is also allowed to vary freely, a lower bound on one scale parameter need not prevent collapse in every direction: a shape matrix could still become nearly singular. The restriction must match the covariance model actually being fitted.

Penalizing Small Covariances

A penalty changes the objective rather than forbidding all covariances below a specified threshold. One direct choice adds a negative inverse-trace term to the log-likelihood:

$$ \ell_{\mathrm{pen}}(\theta) = \ell(\theta) -\frac{\lambda}{2}\sum_{k=1}^K\operatorname{tr}(\Sigma_k^{-1}), \qquad \lambda>0. $$

As a covariance eigenvalue approaches zero, its reciprocal grows without bound and the penalty drives the objective downward. The parameter \(\lambda\) controls the strength of this discouragement; it has units corresponding to a squared measurement scale.

For fixed responsibilities and means, the penalized covariance part for component \(k\) is \[ -\frac{N_k}{2}\log|\Sigma_k| -\frac{1}{2}\operatorname{tr}\{\Sigma_k^{-1}(C_k+\lambda I_d)\}. \] If \(N_k>0\), the same covariance maximization as in the unconstrained Gaussian update gives \(\widehat{\Sigma}_{k,\mathrm{pen}}=(C_k+\lambda I_d)/N_k\). To verify it, set \(B=(C_k+\lambda I_d)/N_k\), which is positive definite. The objective has the form \(-N_k/2\{\log|\Sigma_k|+\operatorname{tr}(\Sigma_k^{-1}B)\}\). Its unique maximum over positive-definite \(\Sigma_k\) is at \(\Sigma_k=B\), as follows by diagonalizing \(B^{-1/2}\Sigma_kB^{-1/2}\) and minimizing \(\log t+1/t\) for each positive eigenvalue \(t\). This proves the update. \(\square\)

Worked Example: An Inverse-Trace Covariance Penalty

Suppose \(N_k=4\) and \(C_k/N_k=\operatorname{diag}(1,0.01)\). Without a penalty, the update is \(\operatorname{diag}(1,0.01)\). Let \(\lambda=0.4\). Then

$$ \widehat{\Sigma}_{k,\mathrm{pen}} = \frac{\operatorname{diag}(4,0.04)+0.4I_2}{4} = \operatorname{diag}(1.1,0.11). $$

Both diagonal entries increase, with the small variance receiving a substantial relative adjustment. For comparison, an eigenvalue floor of \(0.1\) would give \(\operatorname{diag}(1,0.1)\): it leaves the first direction unchanged and sets the second exactly at the threshold. The penalty instead produces an increment of \(\lambda/N_k=0.1\) in each direction for this fixed-responsibility update.

What the Choices Change

MethodWhat is restricted or changed?Typical covariance update
Eigenvalue floorEvery covariance eigenvalue must exceed a fixed threshold.Replace each empirical eigenvalue \(s_j\) by \(\max(s_j,\varepsilon)\).
Structured scale boundA scalar scale in a specified covariance shape must exceed a bound.Replace the unconstrained scale by its maximum with the bound.
Inverse-trace penaltyThe objective is reduced when covariance eigenvalues are small.Add \(\lambda I_d/N_k\) to the empirical covariance in the fixed-responsibility update.

A hard floor has a direct threshold interpretation but can make the fitted solution depend noticeably on the chosen \(\varepsilon\), especially when an empirical eigenvalue lies near it. A structured scale bound is meaningful only alongside its shape assumptions. A penalty gives a continuous adjustment rather than clipping at a boundary, but the amount of adjustment depends on \(\lambda\) and, in the update above, on \(N_k\). These are not interchangeable implementations of the same estimator: each encodes a different restriction or preference.

A positive eigenvalue floor also bounds the height of each Gaussian component density, since \(|\Sigma_k|\geq\varepsilon^d\) implies \[ \phi_d(x;\mu_k,\Sigma_k)\leq (2\pi)^{-d/2}\varepsilon^{-d/2}. \] Consequently the mixture density is bounded by the same constant, and so is its likelihood for a fixed finite data set. This prevents the particular divergence caused by a covariance collapsing to zero. It does not guarantee a unique optimum, good separation of components, or a scientifically useful classification. Likewise, a penalty discourages collapse but does not settle questions about the number or interpretation of components.

Common Pitfall

Wrong reasoning: “A lower bound on component scale is always equivalent to imposing an eigenvalue floor, so either method gives the same estimate.” This overlooks the covariance structure. A scalar scale bound acts on one parameter; an eigenvalue floor acts on every direction of a full covariance matrix. With a fixed shape, the scale bound preserves that shape, while a general eigenvalue floor can raise some directions and leave others unchanged. Moreover, a scale bound alone cannot prevent collapse if an unrestricted shape matrix can still become singular.

A second practical caution is that the displayed EM updates assume \(N_k>0\). An effectively empty component has no ordinary covariance update from its weighted scatter. Implementations need an explicit policy for such components, rather than applying formulas with \(N_k=0\).

Check Your Understanding

Use the definitions and update formulas to answer the following.

  1. If the empirical covariance eigenvalues are \(2\) and \(0.03\), what is the eigenvalue-floored covariance when \(\varepsilon=0.1\), assuming the empirical eigenvectors are unchanged?
  2. Why does an eigenvalue floor prevent a component Gaussian density from becoming arbitrarily high?
  3. For a fixed-shape covariance \(\Sigma=\tau^2A\), what quantity determines the unconstrained estimate of \(\tau^2\)?
  4. In the inverse-trace penalized update, how does the increment to each empirical covariance eigenvalue depend on \(N_k\) and \(\lambda\)?
  5. Why does a scalar scale lower bound fail to guarantee a general eigenvalue floor when the shape matrix is unrestricted?
Clinical Biostats