EM Algorithm seems to work, but Q is not monotonic. Possible reasons?

Question

I have implemented Expectation maximization to fit some of the parameters of a linear Gaussian state space model using Kalman filtering / smoothing.

The model is:

$x(t) = Ax(t - 1) + w(t); w(t) \sim \mathcal{N}(0, Q)$

$y(t) = Hx(t) + v(t); v(t) \sim \mathcal{N}(0, R)$

$m(0) = m_0, P(0|0) = P_0$ // Initial mean and variance

The parameters $\theta$ consist of some of the entries of $A, Q, R$ as well as the initial conditions $m_0, P_0$.

On each iteration of the EM algorithm I (1) set the new initial conditions to $m_0 = \hat{m}(0|T), P_0 = \hat{P}(0 | T)$, where $\hat{m}, \hat{P}$ are the posterior estimates from the smoothing step. and (2) use a numerical optimizer to maximize $\mathcal{Q}(\theta)$ (the expected log likelihood). From what I can see on small 2D examples, this step is correctly computing the maximum of $\mathcal{Q}$. The new initial conditions and new parameters are then fed back into the smoother to do it all again.

To check the correctness of my implementation, I have simulated data from a model with known parameters $\theta_{\text{true}}$ for a variety of sample sizes $N_1, N_2, ...$ (e.g. 100, 1000, ...) and for each of these datasets I run my EM algorithm to convergence to get parameter estimates $\hat{\theta}_1, \hat{\theta}_2, ...$

I am consistently observing that $\hat{\theta}_i \rightarrow \theta_{\text{true}}$, and moreover, this convergence is monotone. So, it seems quite clear that my implementation is correctly inferring the parameters, and it's estimates improve with more data.

I am also looking at the log-likelihood directly, which can be calculated from the filtering recursion, and this function is increasing monotonically. So, running my EM implementation is maximizing the likelihood itself as well.

However, during the EM iterations $t = 1, 2, ...$ the function $\mathcal{Q}(\theta^t)$ is not necessarily monotone in $t$ as it should be. It seems to either increase monotonically (as it should), decrease monotonically, increase and then fall back down before converging, or decrease and then climb back up before converging. In each case however, parameter estimates are still accurate.

Does anyone have suggestions for how this could be the case? I have checked and rechecked my code and my math an innumerable number of times, and am not finding any bugs.

Answer 1

It is the case that the incomplete-data log likelihood has to increase at every step, but is not the case that the expected log likelihood has to increase at every step.

The reason why is hidden in the fact that $\mathcal{Q}(\theta)$ should better be written $\mathcal{Q}(\theta^{t+1}; \theta^t)$. It is true that, as we are choosing the "best" $\theta^{t+1}$ conditional upon the current values $\theta^t$, $\mathcal{Q}(\theta^{t+1}; \theta^t) \geq \mathcal{Q}(\theta^t; \theta^t)$, but it need not be the case that $\mathcal{Q}(\theta^{t+1}; \theta^{t+1}) \geq \mathcal{Q}(\theta^{t+1}; \theta^t)$. As a result:

$$\mathcal{Q}(\theta^{t+1}; \theta^t) \geq \mathcal{Q}(\theta^t; \theta^t) \, ? \, \mathcal{Q}(\theta^t; \theta^{t-1}) \geq \mathcal{Q}(\theta^{t-1}; \theta^{t-1})$$

and no statement can be made about the relationship between $\mathcal{Q}(\theta^{t+1}; \theta^t)$ and $\mathcal{Q}(\theta^t; \theta^{t-1})$.