EM for Gaussian Models with Incomplete Cases

Why Treat Missing Coordinates as Latent Data?

For a multivariate Gaussian model, the likelihood of an incomplete record can be written using only its observed coordinates. Directly optimizing that likelihood, however, may involve several different marginal covariance matrices when records have different missingness patterns. The EM approach instead regards each full record as latent and uses the conditional distribution of its missing coordinates to construct an expected complete-data log likelihood.

The expectation in this construction is taken under the current parameter value, while the candidate parameters in the maximization step remain variable. This distinction is essential: the E-step supplies conditional moments; the M-step maximizes a function built from those moments. An update is defined only if the required maximizer exists.

Setup: Let \(X_1,\ldots,X_n\) be independent \(p\)-dimensional observations with \(X_i\sim N_p(\mu,\Sigma)\), where \(\mu\in\mathbb{R}^p\) and \(\Sigma\) is positive definite. For record \(i\), let \(O_i\) denote its observed coordinates and \(M_i\) its missing coordinates. Assume the missingness mechanism is ignorable for likelihood inference, so the likelihood for \((\mu,\Sigma)\) is the observed-data likelihood. Write the current parameter as \(\theta^{(t)}=(\mu^{(t)},\Sigma^{(t)})\).

Ignorability is an assumption about the missingness mechanism and the likelihood being fitted; it is not a consequence of Gaussianity. If the mechanism is nonignorable, its model generally must be included in the likelihood. The observed-data likelihood uses the Gaussian marginal density for each record's observed coordinates. EM offers a different computational representation of that same likelihood problem.

Complete-Data Log Likelihood and the E-Step

If all coordinates were observed, the Gaussian log likelihood, up to no omitted parameter-dependent terms, would be

\ell_c(\mu,\Sigma;X_1,\ldots,X_n) =-\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma| -\frac{1}{2}\sum_{i=1}^n (X_i-\mu)^\mathsf{T}\Sigma^{-1}(X_i-\mu).

At iteration \(t\), the E-step takes the conditional expectation of this expression given the observed records, using \(\theta^{(t)}\). Define the conditional mean and covariance

m_i^{(t)}=E_{\theta^{(t)}}[X_i\mid X_{i,O_i}=x_{i,O_i}], \qquad V_i^{(t)}=\operatorname{Cov}_{\theta^{(t)}}(X_i\mid X_{i,O_i}=x_{i,O_i}).

The observed coordinates are fixed under this conditional distribution, so their rows and columns in \(V_i^{(t)}\) are zero. If \(O_i\) and \(M_i\) are both nonempty, the conditional Gaussian formula gives

\begin{aligned} E[X_{i,M_i}\mid x_{i,O_i}] &=\mu_{M_i}+\Sigma_{M_iO_i}\Sigma_{O_iO_i}^{-1}(x_{i,O_i}-\mu_{O_i}),\\ \operatorname{Cov}(X_{i,M_i}\mid x_{i,O_i}) &=\Sigma_{M_iM_i} -\Sigma_{M_iO_i}\Sigma_{O_iO_i}^{-1}\Sigma_{O_iM_i}. \end{aligned}

In the E-step these formulas use \(\mu^{(t)}\) and \(\Sigma^{(t)}\). A complete record has \(V_i^{(t)}=0\) and \(m_i^{(t)}=x_i\). A record with no observed coordinates has conditional distribution \(N_p(\mu^{(t)},\Sigma^{(t)})\).

Definition: The EM auxiliary function at iteration \(t\) is \(Q(\mu,\Sigma\mid\theta^{(t)})=E_{\theta^{(t)}}[\ell_c(\mu,\Sigma;X_1,\ldots,X_n)\mid\text{observed data}]\). The expectation is over the missing coordinates, conditional on the observed values. The parameters \((\mu,\Sigma)\) inside \(\ell_c\) are candidate parameters, not the current parameters.

Deriving the Auxiliary Function

For any random vector \(X\) with mean \(m\) and covariance \(V\), and any fixed candidate \(\mu\) and positive-definite \(\Sigma\), the quadratic term satisfies

E[(X-\mu)^\mathsf{T}\Sigma^{-1}(X-\mu)] =\operatorname{tr}(\Sigma^{-1}V)+(m-\mu)^\mathsf{T}\Sigma^{-1}(m-\mu).

To verify this identity, write \(X-\mu=(X-m)+(m-\mu)\), expand the quadratic form, and take expectations. The two cross terms vanish because \(E[X-m]=0\). The remaining random quadratic has expectation \(\operatorname{tr}(\Sigma^{-1}E[(X-m)(X-m)^\mathsf{T}])=\operatorname{tr}(\Sigma^{-1}V)\).

Applying the identity record by record gives the complete expression for the E-step objective:

\begin{aligned} Q(\mu,\Sigma\mid\theta^{(t)}) ={}&-\frac{np}{2}\log(2\pi)-\frac{n}{2}\log|\Sigma|\\ &-\frac{1}{2}\sum_{i=1}^n\left\{ \operatorname{tr}(\Sigma^{-1}V_i^{(t)}) +(m_i^{(t)}-\mu)^\mathsf{T}\Sigma^{-1}(m_i^{(t)}-\mu) \right\}. \end{aligned}

This derivation shows why conditional means alone are not enough: the conditional covariance contributes through the trace term. The M-step maximizes this function over the specified parameter space, holding all \(m_i^{(t)}\) and \(V_i^{(t)}\) fixed. If the maximum exists, a maximizing parameter is the next iterate. If no maximum exists, the stated M-step does not define a finite update.

Worked Example: Conditional Moments for One Partly Observed Record

Suppose the current model is bivariate normal with \(\mu^{(t)}=(0,0)^\mathsf{T}\) and \(\Sigma^{(t)}=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}\). The first coordinate is observed as \(x_1=2\), while the second is missing. Find the conditional moments used in the E-step.

The observed-coordinate variance is \(4\), the covariance between coordinates is \(1.2\), and the missing-coordinate variance is \(1\). Thus the conditional mean of coordinate 2 is \(0+(1.2/4)(2-0)=0.6\).

Its conditional variance is \(1-(1.2)(1/4)(1.2)=1-0.36=0.64\). The first coordinate remains fixed at 2, so the full conditional mean and covariance are \(m^{(t)}=(2,0.6)^\mathsf{T}\) and \(V^{(t)}=\begin{pmatrix}0&0\\0&0.64\end{pmatrix}\). The zero in the first diagonal entry records that the observed value has no conditional uncertainty.

Worked Example: Building Expected Sufficient Statistics

Let \(X_i\sim N_2(0,\Sigma^{(t)})\), where \(\Sigma^{(t)}=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\). There are two records: the first is fully observed as \((1,0)^\mathsf{T}\), and the second has only its first coordinate observed, at \(-1\). Compute the sum of the conditional expectations \(E[X_iX_i^\mathsf{T}\mid\text{observed data}]\).

For the complete record, \(m_1=(1,0)^\mathsf{T}\), \(V_1=0\), and \(E[X_1X_1^\mathsf{T}\mid X_1]=(1,0)^\mathsf{T}(1,0)=\begin{pmatrix}1&0\\0&0\end{pmatrix}\).

For the second record, the conditional mean of coordinate 2 is \(0.5(-1)=-0.5\), and its conditional variance is \(1-0.5^2=0.75\). Hence \(m_2=(-1,-0.5)^\mathsf{T}\) and \(V_2=\begin{pmatrix}0&0\\0&0.75\end{pmatrix}\). Since \(E[XX^\mathsf{T}]=V+mm^\mathsf{T}\) for any vector with mean \(m\) and covariance \(V\), \(E[X_2X_2^\mathsf{T}\mid X_{2,1}=-1]=\begin{pmatrix}1&0.5\\0.5&0.25+0.75\end{pmatrix}=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\).

Adding the two contributions gives \(\sum_i E[X_iX_i^\mathsf{T}\mid\text{observed data}]=\begin{pmatrix}2&0.5\\0.5&1\end{pmatrix}\). This is a conditional expectation under the current model, not an assertion that either missing value was observed.

Worked Example: Why Single Imputation Changes the E-Step

Use the second record of the preceding example, with conditional mean \((-1,-0.5)^\mathsf{T}\) and conditional variance \(0.75\) for the missing coordinate. Evaluate its expected quadratic contribution to \(Q\) at candidate \(\mu=(0,0)^\mathsf{T}\) and \(\Sigma=I_2\).

The correct expected quadratic is \(\operatorname{tr}(I_2V)+m^\mathsf{T}m=0.75+[(-1)^2+(-0.5)^2]=0.75+1.25=2\).

If instead the missing coordinate is replaced by its conditional mean and treated as known, the quadratic would be \(1.25\). The difference, \(0.75\), is precisely the conditional variance contribution. Because the complete-data log likelihood includes minus one half of this quadratic, the correct contribution is \(-1\), rather than the plug-in value \(-0.625\). A conditional mean imputation is therefore not a substitute for taking the conditional expectation of the log likelihood.

The EM Iteration and Its Existence Condition

For a general parameter space, one iteration consists of computing the conditional moments under \(\theta^{(t)}\), forming \(Q(\theta\mid\theta^{(t)})\), and choosing a maximizer over the allowed candidate parameters if such a maximizer exists. In a constrained or structured Gaussian model, the maximization is over that model's permitted parameter set; one must check existence there rather than assume it.

A useful exact check is available for the unrestricted Gaussian model. Define the conditional expected centered scatter at a candidate mean \(\mu\) by

A^{(t)}(\mu)=\sum_{i=1}^n\left\{V_i^{(t)}+ (m_i^{(t)}-\mu)(m_i^{(t)}-\mu)^\mathsf{T}\right\}.

This matrix is positive semidefinite, and the parameter-dependent part of \(Q\) is \(-\frac n2\log|\Sigma|-\frac12\operatorname{tr}(\Sigma^{-1}A^{(t)}(\mu))\). The unrestricted M-step has a finite maximizer precisely when the minimized expected centered scatter is positive definite. Here is the reason.

Center at the conditional means.
Let \(\bar m^{(t)}=n^{-1}\sum_i m_i^{(t)}\). Expanding around this value gives \(A^{(t)}(\mu)=A^{(t)}(\bar m^{(t)})+n(\mu-\bar m^{(t)})(\mu-\bar m^{(t)})^\mathsf{T}\). The added matrix is positive semidefinite, so the smallest scatter in the positive-semidefinite ordering occurs at \(\bar m^{(t)}\). Write \(C=A^{(t)}(\bar m^{(t)})\).

Check positive definiteness.
If \(C\) is positive definite, set \(B=C/n\) and \(T=B^{1/2}\Sigma^{-1}B^{1/2}\). Apart from a constant, the objective at \(\mu=\bar m^{(t)}\) is \(\frac n2\{\log|T|-\operatorname{tr}(T)\}\). For each positive eigenvalue \(\lambda\) of \(T\), \(\log\lambda-\lambda\leq -1\), with equality only at \(\lambda=1\). Thus a finite maximum is attained at \(T=I_p\), equivalently \(\Sigma=B\).

Recognize the singular case.
If \(C\) is singular, choose an orthonormal basis in which one coordinate is in its null space, and let the corresponding covariance eigenvalue tend to zero while keeping the others fixed and positive. The trace term involving \(C\) stays finite, while \(-\frac n2\log|\Sigma|\) tends to \(+\infty\). Therefore the unrestricted auxiliary function is unbounded above: there is no finite M-step maximizer.

This argument is an existence check, not a guarantee that every other model's M-step has a solution. For the unrestricted Gaussian M-step, the required condition is that \(C\), the expected centered scatter at its minimizing center, be positive definite. If it fails, suitable constraints or regularization may be needed to define a finite update.

Common Pitfalls and Practical Meaning

The latent-record construction is valuable because it expresses each E-step through conditional moments and leaves the M-step as an optimization of a complete-data-form objective. It also makes clear which uncertainty is being integrated over. With different missingness patterns, each record has its own conditional distribution, but all contributions enter the same \(Q\) function.

Plugging in conditional means. This drops \(V_i^{(t)}\) and changes the expected complete-data log likelihood. Include both conditional means and covariances.
Updating the moments during maximization. The E-step moments are computed using \(\theta^{(t)}\) and held fixed throughout the M-step. Recomputing them as candidate parameters change is not the defined EM auxiliary-function maximization.
Assuming every M-step exists. An unconstrained Gaussian auxiliary function can be unbounded above when its expected centered scatter is singular. Do not label an arbitrary parameter value as a maximizer in that case.
Ignoring the missingness assumption. The Gaussian-data EM construction above uses an ignorable missingness mechanism. A nonignorable mechanism can require additional latent variables and likelihood terms.

Key takeaway: For incomplete Gaussian records, the E-step replaces each full-data quadratic contribution by its conditional expectation under the current parameters. The resulting \(Q\) function contains both conditional means and conditional covariances; an EM update requires an actual maximizer over the chosen parameter space.

Check Your Understanding

Use the definitions and derivations above to answer the following.

In the E-step, which parameters determine the conditional distribution of the missing coordinates, and which parameters are varied in the M-step?
Why does the conditional covariance appear in the expected Gaussian quadratic term?
For a complete record, what are its conditional mean and covariance given the observed data?
What condition on the unrestricted Gaussian expected centered scatter ensures a finite M-step maximum?
Explain why replacing a missing value by its conditional mean does not generally compute the EM auxiliary function.