# What is the M-step for Expectation Maximization for a multivariate Gaussian hidden Markov model with missing observations?

Is there a closed-form solution for the M-step of the expectation maximization (EM) algorithm for a multivariate Gaussian hidden Markov model (HMM) where observations are missing at random? $\newcommand\given[1][]{\:#1\vert\:}$

My HMM notation, largely borrowed from https://www.princeton.edu/~rvan/orf557/hmm080728.pdf, is:

• $X_t \in \mathrm{E}$ is a discrete hidden state with transition matrix $P$ and initial distribution $\mu$
• The state space $\mathrm{E} \equiv \left\{1, 2, \ldots, d\right\}$ is finite
• The observations are $\left(\,Y_t \given X_t = i\,\right)\sim \mathcal{N}\left(\mathbf{m}_i, \mathbf{v}_i\right)$
• If the observations are scalars, $\mathbf{m}$ and $\mathbf{v}$ are simply vectors of means and variances (one for each hidden state)
• If the observations are multivariate normal, say $Y_t \in \mathbb{R}^k$, each $\mathbf{m}_i$ is a $k$-vector of means and each $\mathbf{v}_i$ is a $k$-by-$k$ covariance matrix

I'm interested in estimating $P$, $\mu$, $\mathbf{m}$ and $\mathbf{v}$ using the EM algorithm. Let $\mathbf{y}$ denote an $n$-by-$k$ matrix of observed data. Let $\pi_{it}$ denote the posterior probability over hidden state $i$ for observation $t$, computed under the previous iteration's parameters, i.e. $\pi_{it} = \Pr[X_t = i \given \mathbf{y}]$.

The notes I linked to earlier give the M step when $k$ = 1, i.e. scalar Gaussian observations (section 6.2, page 87):

$$\begin{equation*} \mathbf{m}_i = \frac{\sum^n_{t=1} \mathbf{y}_t \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$ $$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right)^2 \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$

I'm interested mainly in the updates to $\mathbf{m}$ and $\mathbf{v}$ so I'll omit the updates to $P$ and $\mu$ for clarity.

The M-step updates for multivariate Gaussian observations are quite similar. The update to $\mathbf{m}_i$ is identical (the only difference being that $y_t$ is now a $k$-vector), and the update for the covariance matrices is:

$$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right) \, \left(\mathbf{y}_t - \mathbf{m}_i\right)^{\intercal} \, \pi_{it}}{\sum^n_{t=1} \pi_{it}} \end{equation*}$$

This equation appears on slide 9 of http://personal.ee.surrey.ac.uk/Personal/P.Jackson/tutorial/hmm_tut4.pdf, for example, with slightly different notation.

What is the M-step if elements of $\mathbf{y}$ (the $n$-by-$k$ matrix of observations) are missing at random? My guess is to first define a matrix $\tilde{\mathbf{y}}$ containing 1s where $\mathbf{y}$ is observed and 0s where it is missing (i.e. a matrix of indicators for non-missing data). I then recode the missing observations in $\mathbf{y}$ as 0.

My updates (after the recoding described above) are:

$$\begin{equation*} \mathbf{m}_i = \frac{\sum^n_{t=1} \mathbf{y}_t \, \pi_{it}}{\sum^n_{t=1} \tilde{\mathbf{y}}_t \, \pi_{it}} \end{equation*}$$ $$\begin{equation*} \mathbf{v}_i = \frac{\sum^n_{t=1} \left(\mathbf{y}_t - \mathbf{m}_i\right) \, \left(\mathbf{y}_t - \mathbf{m}_i\right)^{\intercal} \, \pi_{it}}{\sum^n_{t=1} \tilde{\mathbf{y}}_t \, \tilde{\mathbf{y}}_t^{\intercal} \pi_{it}} \end{equation*}$$

In the equation above, both the numerator and denominator are $k$-by-$k$ matrices. The fraction should be read as an element-wise division.

Is this correct -- is it even guaranteed to produce a positive definite covariance matrix? Does a correct analytical update exist? Any pointers to papers or lecture notes that discuss multivariate Gaussian HMMs with observations missing at random would be much appreciated.

A relevant paragraph from https://en.wikipedia.org/wiki/Estimation_of_covariance_matrices:

In the general case, the unbiased estimate of the covariance matrix provides an acceptable estimate when the data vectors in the observed data set are all complete: that is they contain no missing elements. One approach to estimating the covariance matrix is to treat the estimation of each variance or pairwise covariance separately, and to use all the observations for which both variables have valid values. Assuming the missing data are missing at random this results in an estimate for the covariance matrix which is unbiased. However, for many applications this estimate may not be acceptable because the estimated covariance matrix is not guaranteed to be positive semi-definite. This could lead to estimated correlations having absolute values which are greater than one, and/or a non-invertible covariance matrix.

Do any references discuss how to deal with this problem in the context of multivariate Gaussian HMMs?

Here's a little R example showing covariance matrix estimates with negative eigenvalues:

library(mixtools)  # For rmvnorm

get_sigma_hat <- function(n_obs, sigma, pr_missing) {
## Generate random multivariate data y with observations missing at random
## Return estimated variance-covariance matrix
stopifnot(0 <= pr_missing && pr_missing < 1)
stopifnot(all(eigen(sigma)$values > 0.0)) stopifnot(identical(sigma, t(sigma))) y <- rmvnorm(n_obs, sigma=sigma) stopifnot(all(dim(y) == c(n_obs, nrow(sigma)))) y_non_missing <- matrix(1*(runif(n_obs * nrow(sigma)) > pr_missing), n_obs, nrow(sigma)) mean(y_non_missing) # Should be close to (1 - pr_missing) observed_y <- ifelse(y_non_missing, y, 0) # Careful, NAs recoded to zero sigma_hat <- t(observed_y) %*% observed_y / (t(y_non_missing) %*% y_non_missing) return(list(sigma_hat=sigma_hat, sigma_hat_positive_definite=all(eigen(sigma_hat)$values > 0.0)))
}

set.seed(543987)

sigma <- matrix(c(1.0, 0.5, 0.6,
0.5, 2.0, 0.7,
0.6, 0.7, 3.0), 3, 3)  # True variance-covariance matrix for observations in R^3
n_reps <- 1000
reps <- replicate(n_reps, get_sigma_hat(n_obs=500, sigma=sigma, pr_missing=0.10), simplify=FALSE)

mean(sapply(reps, function(x) x$sigma_hat_positive_definite)) # 1.0 i.e. all positive definite mean(sapply(reps, function(x) x$sigma_hat[1, 1]))  # Close to sigma[1, 1]
mean(sapply(reps, function(x) x$sigma_hat[3, 2])) # Close to sigma[3, 2] reps <- replicate(n_reps, get_sigma_hat(n_obs=400, sigma=sigma, pr_missing=0.80), simplify=FALSE) mean(sapply(reps, function(x) x$sigma_hat_positive_definite))  #  Some negative eigenvalues!
sigma_hat <- Filter(function(x) !x$sigma_hat_positive_definite, reps)[[1]]$sigma_hat
eigen(sigma_hat)$values # Third eigenvalue is negative cov2cor(sigma_hat) # ...but all correlations are in [-1, 1] another_sigma_hat <- Filter(function(x) !x$sigma_hat_positive_definite, reps)[[2]]$sigma_hat eigen(another_sigma_hat)$values  # Again, third eigenvalue is negative
cov2cor(another_sigma_hat)  # This one has some correlations above 1


The code simulates a simplified version of the problem I describe above, where the $Y$s have mean zero and there are no weights $\pi$.