Stable Log-Likelihood and Responsibility Computation

Why Mixture Calculations Need Numerical Care

A Gaussian mixture assigns each observation a density that is a weighted sum of component densities. In moderate dimensions, or for observations far from all component means, those densities can be extremely small. A computer may round them to zero even though their logarithms are representable. Taking the logarithm after that rounding produces \(-\infty\), losing information about the likelihood and potentially making posterior component probabilities undefined.

The remedy is to do the calculation on the log scale. The key operation is the log-sum-exp: a stable way to compute the logarithm of a sum of exponentials without first forming tiny exponentials. The same calculation gives stable posterior responsibilities.

For computation, define the component log contribution \(a_{ik}=\log\pi_k+\log f_k(x_i)\), with \(\log 0=-\infty\). Then the contribution of observation \(i\) to the log-likelihood is \(\log\sum_k\exp(a_{ik})\). The responsibilities are obtained by normalizing those same component contributions.

The Shifted Log-Sum-Exp Calculation

Exponentiating several large negative numbers is the source of underflow. Instead, shift all values by their maximum. For finite \(a_1,\ldots,a_K\), put \(m=\max_k a_k\). Since \(a_k-m\leq 0\), the exponentials \(\exp(a_k-m)\) are at most one, and at least one equals one.

Proof. Factoring \(e^m\) out of the sum gives \(\sum_k e^{a_k}=e^m\sum_k e^{a_k-m}\). Taking logarithms gives the stated identity. Also, \(e^{a_k-s}=e^{a_k}/\sum_j e^{a_j}\). Summing this expression over \(k\) gives one, and each term is nonnegative. \(\square\)

For a mixture observation, calculate \(m_i=\max_k a_{ik}\), then use $$ \log f(x_i)=m_i+\log\left(\sum_{k=1}^K e^{a_{ik}-m_i}\right),\qquad \tau_{ik}=\exp\left(a_{ik}-m_i-\log\left(\sum_{j=1}^K e^{a_{ij}-m_i}\right)\right). $$

These formulas do not change the mathematical answer; they change the order of arithmetic so that the computer can retain it. Responsibilities can be very small, but their exponent arguments are never positive. If some component has zero weight, its log contribution is \(-\infty\) and its responsibility is zero. If all contributions for an observation are \(-\infty\), then the mixture density there is zero: the log-likelihood contribution is genuinely \(-\infty\), and responsibilities are not defined by the usual ratio. That case is different from numerical underflow.

Computing Gaussian Log Densities

Log-sum-exp stabilizes the mixture sum, but each component log density must also be calculated reliably. For \(X\sim N(\mu,\Sigma)\) in \(d\) dimensions, assume \(\Sigma\) is symmetric positive definite. The log density at \(x\) is $$ \log f(x)=-\frac d2\log(2\pi)-\frac12\log|\Sigma|-\frac12(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu). $$

A numerical implementation should generally avoid explicitly forming \(\Sigma^{-1}\). Let \(\Sigma=LL^\mathsf{T}\) be its Cholesky factorization, with \(L\) lower triangular, and solve \(Ly=x-\mu\). Then the quadratic term is \(y^\mathsf{T}y\), and \(\log|\Sigma|=2\sum_{j=1}^d\log L_{jj}\). Thus the calculation uses a triangular solve and logarithms of positive diagonal entries rather than a matrix inverse and determinant.

Responsibilities and the Common Shift

Responsibilities are posterior probabilities, not merely a convenient normalization. Suppose \(\pi_k f_k(x_i)\) are the joint density contributions for \(X_i=x_i\) and \(Z_i=k\). The marginal density of \(X_i\) is their sum. Bayes’ rule therefore gives the responsibility as the component contribution divided by that sum. Substituting \(a_{ik}\) yields the stable expression above.

Proof. By definition, \(e^{a_{ik}}=\pi_k f_k(x_i)\). Let \(s_i=\log\sum_j e^{a_{ij}}\), so \(e^{s_i}=\sum_j\pi_j f_j(x_i)=f(x_i)>0\). Hence \(e^{a_{ik}-s_i}=\pi_k f_k(x_i)/f(x_i)\), which is Bayes’ rule for \(P(Z_i=k\mid X_i=x_i)\). Summing over \(k\) gives \(\sum_k \tau_{ik}=\sum_k\pi_k f_k(x_i)/f(x_i)=1\). \(\square\)

Applying the Calculation to a Data Set

For each observation, form its \(K\) log contributions, compute a shifted log-sum-exp, and add the resulting log contributions to obtain \(\ell\). A numerically stable implementation should never first compute the component densities and then take their logarithms. It should compute component log densities, include log weights, and normalize the component contributions on the log scale.

Common Mistakes and Practical Checks

A frequent mistake is to calculate \(\sum_k\pi_kf_k(x_i)\) directly, then take its logarithm. If all component densities round to zero, the result is \(-\infty\), even when the true log-likelihood is finite. The correction is to calculate \(a_{ik}\) first and use the shifted log-sum-exp formula.

Another mistake is to exponentiate \(a_{ik}\) to obtain responsibilities and divide by their direct sum. When every term underflows, this becomes \(0/0\). Instead compute each \(\exp(a_{ik}-m_i)\), sum those shifted values, and divide by their sum—or equivalently use the log responsibility formula.

Finally, a stable summation cannot repair an invalid component density calculation. Check that covariance matrices meet the assumptions of the Gaussian density formula, and use a stable factorization rather than relying on a numerically fragile inverse. Useful checks are that every finite observation-level mixture log-density is finite when all component Gaussian densities are positive, responsibilities are nonnegative, and their sum is one up to rounding error.