How to use log probabilities in PCA mixture EM algorithm

I'm trying to implement PCA mixtures (Tipping & Bishop 2006 Appendix C) on the Tobomovirus. I'll summarize the mathematical background and algorithm here:

For a single PCA model, we assume a latent variable model

$$\mathbf{t} = \mathbf{Wx} + \boldsymbol\mu + \boldsymbol\epsilon.$$ This implies a probability distribution $$p(\mathbf{t}|\mathbf{x}) = (2\pi\sigma^2)^{-d/2}\mathrm{exp}\{-\frac{1}{2\sigma^2}||\mathbf{t} - \mathbf{Wx}- \boldsymbol\mu||^2\}$$ since the noise term $$\boldsymbol\epsilon$$ is zero-mean Gaussian with variance $$\sigma^2$$.

We define the Gaussian prior

$$p(\mathbf{x}) = (2\pi)^{-q/2}\mathrm{exp}\{-\frac{1}{2}\mathbf{x}^T\mathbf{x}\}$$ which allows us to derive the marginal distribution

$$p(\mathbf{t}) = (2\pi)^{-d/2}|\mathbf{C}|^{-1/2}\mathrm{exp}\{-\frac{1}{2}(\mathbf{t} - \boldsymbol\mu)^T\mathbf{C}^-1(\mathbf{t} - \boldsymbol\mu)\}$$

where $$\mathbf{C} = \sigma^2\mathbf{I} + \mathbf{WW}^T$$.

The PCA mixture model means that we use a mixture of PCA models like this:

$$p(\mathbf{t}) = \sum_{i=1}^M \pi_i p(\mathbf{t}|i)$$

The algorithm:

1. Initialize $$\pi_i$$, $$\boldsymbol\mu_i$$, $$\mathbf{W}_i$$, $$\sigma_i^2$$. These are the mixing coefficient, mean (of the marginal), weight matrix and variance of each model.
2. Calculate the responsibilities,

$$R_{ni} = \frac{p(\mathbf{t}_n|i)\pi_i}{p(\mathbf{t}_n)}$$ where $$p(\mathbf{t}|i)$$ is the probability for a single PCA model $$i$$.

1. Calculate mixing coefficients

$$\tilde\pi_i = \frac{1}{N}\sum^N_{n=1} R_{ni}.$$

1. Calculate the means

$$\tilde{\boldsymbol\mu}_i = \frac{\sum_{n=1}^N R_{ni}\mathbf{t}_n}{\sum_{n=1}^NR_{ni}}$$

1. Calculate the new weight matrices

$$\mathbf{\tilde W}_i = \mathbf{S}_i \mathbf{W}_i(\sigma^2 \mathbf{I}+\mathbf{M}^{-1}\mathbf{W}_i^T\mathbf{S}_i\mathbf{W}_i)^{-1}$$

where $$\mathbf{S}_i = \frac{1}{\tilde\pi_i N}\sum^N_{n=1}R_{ni}(\mathbf{t}_n - \tilde{\boldsymbol\mu})(\mathbf{t}_n - \tilde{\boldsymbol\mu})^T$$.

1. Repeat from step 2.

The tilde operator means that it is a new parameter being calculated.

Here's my problem: When implementing this on the Tobomovirus dataset, the long 18x1 feature vectors and the resulting 18x18 $$\mathbf{C}$$ matrix lead to the probabilities $$p(\mathbf{t})$$ in step 2 becoming very very small, leading to underflow. I would like to work with the log probabilities instead, but I'm not sure how this will affect the algorithm and what exactly I need to change.

I know nothing about "PCA mixtures EM algorithms", but I can answer how I solve this problem with general EM algorithms, such as mixture of multivariate Bernoulli models. Usually, the EM algo involves calculating terms of the form

$$P(z|\underline{x}, \underline{\theta})$$ in which z is your latent variable, $$\underline{x}$$ are your features and $$\underline{\theta}$$ are your parameters which you're trying to optimise your log-likelihood wrt. This is also usually the most challenging stage of the algorithm, from a numerical standpoint.

Assuming your latent variable is discrete:

$$P(z|\underline{x}, \underline{\theta})=\frac{P(\underline{x}|z, \underline{\theta})}{\sum_{k}P(\underline{x}| z=k, \underline{\theta})}$$

Numerical instabilities occur because if you have many features, $$P(\underline{x}|z, \underline{\theta})$$ will always be small, so you are dividing a small number by the sum of many other small numbers.

One way to mitigate this is to calculate $$P^{*}=max_{k}\left[P(\underline{x}|z=k, \underline{\theta})\right]$$ , and then write

$$P(z|\underline{x}, \underline{\theta})=\frac{P(\underline{x}|z, \underline{\theta})}{P^{*}}\left(\sum_{k}\frac{P(\underline{x}|z=k, \underline{\theta})}{P^{*}}\right)^{-1}$$

These ratios are easier to calculate if you have a way of calculating $$\log P(\underline{x}|z, \underline{\theta})$$ which does not involve directly evaluating $$P(\underline{x}|z, \underline{\theta})$$. For example, if your likelihood function stipulates that all features are independent on one and other, than the likelihood of any data point is a product over features, thus the log likelihood is a sum over logarithms and that will be easy to evaluate with no numerical instabilities. Likewise if your likelihood is some exponential family, taking the logarithm is the same as evaluating the exponent which should be nice and finite.

To calculate $$\frac{P(\underline{x}|z, \underline{\theta})}{P^{*}}$$, you can evaluate $$e^{\log P(\underline{x}|z, \underline{\theta}) - \log P^{*}}$$ instead

If there is a big difference between $$max_{k}\left[P(\underline{x}|z=k, \underline{\theta})\right]$$ and $$min_{k}\left[P(\underline{x}|z=k, \underline{\theta})\right]$$, and it can happen that the difference is more orders of magnitude than the difference between the largest and smallest float in python, in my experience, then for some j, $$\frac{P(\underline{x}|z=j, \underline{\theta})}{P^{*}}$$ will be zero, even when using the logarithm trick, which is fine, it's the algorithm telling you that there is effectively zero probability that cluster j (or latent variable value j) created this data point.

Is this useful for a "PCA mixture EM algorithm" ?

• Hi, thanks for your answer. How can I calculate $P(z=k|x,θ)$ to maximize $P^*$ when not being able to calculate $P(z=k|x,θ)$ is my main problem? Did you perhaps mean to write $P^*=max_k[P(x|z=k,θ)]$? – Sahand Jan 10 at 12:59
• there was a typo in the original. I've made some edits. Does the expression for $P^{*}$ make sense now? – gazza89 Jan 10 at 13:52
• Yes indeed! This might in fact solve the problem, I'm implementing it now. Thanks again – Sahand Jan 10 at 13:56
• It works perfectly. – Sahand Jan 11 at 17:37