# MAP estimation for a Gaussian mixture using EM. Concerns with the covariance update formula

I am implementing the EM algorithm for a Gaussian mixture model with prior; that is, I am using the EM algorithm to find the MAP estimate, rather than the ML estimate. As briefly discussed in section 9.4 of Bishop's Pattern Recognition and Machine Learning this can be done by modifying the M-step to optimize

$$R(\theta, \theta^\text{old}) = Q(\theta, \theta^\text{old}) + \ln p(\theta) \tag 1$$

instead of $$Q(\mathbf \theta, \mathbf \theta^\text{old})$$, where $$Q(\theta, \theta^\text{old})$$ is familiar from the regular EM algorithm. Now, I think I have correctly worked out the parameter values that maximize (1), but I am a bit concerned about the result for the component covariance matrices.

Recall that a Gaussian mixture distribution has the form

$$p(x | \pi, \mu, \Lambda) = \sum_{k=1}^K \pi_k \mathcal N(x | \mu_k, \Lambda_k^{-1}), \tag 2$$

with the mixing coefficients $$\pi = \{\pi_k\}$$, component means $$\mu = \{\mu_k\}$$, and component precision matrices $$\Lambda = \{\Lambda_k\}$$ being the parameters of the distribution. We treat these parameters as latent variables, and try to determine them through EM on a data set $$X = \{x_n\}$$.

The conjugate prior for the mixing components is a Dirichlet distribution,

$$p(\pi) = \text{Dir}(\pi|\alpha) = C(\alpha) \prod_{k=1}^K \pi_k^{\alpha_k - 1}, \tag 3$$

and the conjugate prior for each mean-precision pair is a Gaussian-Wishart distribution, giving

$$p(\mu, \Lambda) = \prod_{k=1}^K p(\mu_k, \Lambda_k) = \prod_{k=1}^K \mathcal N(\mu_k | \mu_0, (\beta_0 \Lambda_k)^{-1}) \mathcal W(\Lambda_k | W_0, \nu_0). \tag 4$$

Optimizing (1) with respect to $$\mu$$ or $$\Lambda$$ is simply done by setting derivatives to zero and solving. To optimize with respect to $$\pi$$ one must employ a Lagrange multiplier to enforce the condition $$\sum_{k=1}^K \pi_k = 1$$. Assuming I have not made any mistakes, the result is as follows:

$$\pi_k^* = \frac{N_k + \alpha_k - 1}{N + \sum_{j=1}^K \alpha_j - K} \tag 5$$

$$\mu_k^* = \frac{\sum_{n=1}^N \gamma_{nk} x_n + \beta_0 \mu_0}{N_k + \beta_0} \tag 6$$

$$\Sigma_k^* = \frac{1}{N_k + \nu_0 - D} \left[ \sum_{n=1}^N \gamma_{nk} (x_n - \mu_k^*)(x_n - \mu_k^*)^T \\+ \beta_0 (\mu_k^* - \mu_0)(\mu_k^* - \mu_0)^T + W_0^{-1} \right] \tag 7$$

Here $$\{\gamma_{nk}\}$$ are the so called responsibilities

$$\gamma_{nk} = \frac{\pi_k \mathcal N(x_n|\mu_k,\Lambda_k^{-1})}{\sum_{j=1}^K \pi_j \mathcal N(x_n|\mu_j,\Lambda_j^{-1})} \tag 8$$

which were evaluated in the E step (with respect to the old parameter values), and $$N_k = \sum_{n=1}^N \gamma_{nk}$$ (sometimes viewed as an effective count of points assigned to component $$k$$). $$D$$ is the dimensionality of the data points $$x_n$$.

So, finally to my question. It seems to me that the denominator $$N_k + \nu_0 - D$$ in (7) may in some cases become negative, which would imply that $$\Sigma_k^*$$ fails to be positive definite, as is required for a covariance matrix. Say for example that we initialize the $$k$$th component such that $$\mu_k$$ is very far away from all of the data. Then $$\gamma_{nk}$$ would be vanishingly small for all $$n$$, and hence $$N_k \ll 1$$. From the definition of the Wishart distribution we have $$\nu_0 > D-1$$, but this is not enough to guarantee $$N_k + \nu_0 - D > 0$$, since $$\nu_0$$ is a continuous parameter that can be arbitrarily close to $$D-1$$.

How should I interpret this result? Does a maximum of (1) not exist in this sort of case? Have I simply made a mistake in my calculation? Is there perhaps a written source somewhere that I can check my results against?

• Practically, you could bound the hyperparameters Oct 7, 2023 at 3:19
• Since $\nu_0$ is a parameter of the prior,ie an hyperparameter, it can be chosen such that $\nu_0\ge D$. This is a reasonable constraint when considering that, when $\nu_0<D$, the Whishat matrix $\Lambda_k$ is not invertible with probability $1$.. Oct 7, 2023 at 7:58
• Furthermore, the existence of a $-k$-th order moment is constrained by$$\nu_0>2k+D-1$$which means that there is no (prior) expectation of the Wishart matrix when $\nu_0\le D+1$. Oct 7, 2023 at 8:09
• @Xi'an Those are excellent points, both of which I did not know. Still, I am hesitant to constrain $\nu_0$ further in my implementation, since the Wishart distribution is well defined for $\nu_0 > D - 1$, despite the issues you raise. Though in practice I would heed your advice when actually using the implementation.
– ummg
Oct 7, 2023 at 16:08
• @Xi'an Feel free to write an answer if you want to get credited for your comments, which have been very helpful, and I will mark it as accepted. I might also write an answer of my own.
– ummg
Oct 7, 2023 at 19:39

While this is not a purely Bayesian framework (since MAP estimators can be argued as being un-Bayesian), one should consider that the prior $$p(\theta)$$ in this analysis is a genuine prior, thus chosen by a Bayesian experimenter rather than provided by Nature. This freedom of choice (and the ensuing subjectivity of the Bayesian analysis) mean that constraints to ensure the convergence of the modified EM sequence can be enforced upon the prior by picking the hyperparameter $$\nu_0$$ accordingly, i.e. $$\nu_0>D.$$
Note that there are several arguments for choosing $$\nu_0$$ larger than the lower bound $$D-1$$:
• when $$𝜈_0<𝐷$$, the Whishart random matrix $$Λ_𝑘$$ is not almost surely invertible, i.e., the covariance matrix may be degenerate
• when $$𝜈_0\ge𝐷$$, the Whishart random matrix $$Λ_𝑘$$ is thus almost surely positive definite
• when $$𝜈_0>2n+𝐷−1$$ the inverse Whishart random matrix $$Λ_𝑘^{-1}$$ enjoys an $$n$$-th order moment, which implies that
• when $$𝜈_0≤𝐷+1$$ there is no (prior) expectation of the inverse Wishart random matrix $$Λ_𝑘^{-1}$$.
As pointed out in a comment (and later an answer) by Xi'an the matrix $$\Lambda_k \sim \mathcal W(W_0, \nu_0)$$ is invertible (and hence positive definite) with probability 1 if $$\nu_0 \ge D$$, but not otherwise. Hence, when used as a prior for a precision matrix, it makes little sense to let $$\nu_0 < D$$ in the Wishart distribution.
Moreover, it is not possible to salvage (7) for the case $$N_k + \nu_0 - D < 0$$ by performing a constrained optimization that enforces the positive definiteness of $$\Lambda_k$$ (or $$\Sigma_k$$). To see this, note that $$R(\theta, \theta^\text{old})$$ contains a term of the form $$(N_k + \nu_0 - D) \ln |\Lambda_k|$$. This term can be made arbitrary large by making $$\Lambda_k$$ sufficiently close to singular. Hence $$R(\theta, \theta^\text{old})$$ cannot have a maximum on the set of positive definite $$\Lambda_k$$.