EM Algorithm for Gaussian Factor Analysis

Why Use EM for Gaussian Factor Analysis?

In a Gaussian factor model, the observed variables are related through latent factors, while each observed variable also has its own specific error. The likelihood is therefore built from a covariance matrix of the form \(\Lambda\Lambda^\mathsf{T}+\Psi\). Direct likelihood optimization can be difficult because the factors are unobserved. The expectation-maximization (EM) algorithm handles this by alternating between estimating the conditional distribution of the factors and maximizing an expected complete-data log-likelihood.

The latent-factor representation makes both steps explicit. In the E-step, the current parameters give the conditional means and covariances of the factors given the observations. In the M-step, those conditional moments act like sufficient statistics for updating the loadings and uniquenesses. The mean deserves particular care: if it is unknown, its EM update is generally not just the observed sample mean.

The E-Step: Conditional Factor Moments

Write the current parameter values as \(\mu,\Lambda,\Psi\), and let \(x_i\) denote the observed value of \(X_i\). Completing the square in the joint Gaussian density of \(Z_i\) and \(X_i\) gives the conditional distribution of the factor:

The E-step needs the first and second conditional moments. They are \(E[Z_i\mid x_i]=m_i\) and \(E[Z_iZ_i^\mathsf{T}\mid x_i]=V+m_im_i^\mathsf{T}\). The second identity follows by writing \(Z_i=(Z_i-m_i)+m_i\): the centered conditional term has covariance \(V\), and its conditional mean is zero.

The M-Step: Loadings, Mean, and Uniquenesses

Let \(m_i\) and \(V\) be the E-step quantities computed at the current parameters. Define the average observed vector and average conditional factor mean by \(\bar x=n^{-1}\sum_i x_i\) and \(\bar m=n^{-1}\sum_i m_i\). The following centered expected moments are useful:

Here \(A\) is \(q\times q\) and \(B\) is \(p\times q\). In the ordinary unconstrained M-step, the updates are \(\Lambda_{\mathrm{new}}=BA^{-1}\) and \(\mu_{\mathrm{new}}=\bar x-\Lambda_{\mathrm{new}}\bar m\). Since \(V\) is positive definite under the stated assumptions, \(A\) is positive definite and its inverse exists.

To derive the loading and mean updates, retain the terms in the expected complete-data log-likelihood that depend on them. Apart from constants, these terms are minus one half the sum of the expected squared residuals weighted by \(\Psi^{-1}\). Because \(\Psi\) is diagonal, each observed-variable row minimizes its own expected squared residual. For a proposed \(\mu\), the rowwise least-squares normal equations give \(\Lambda=\{\sum_i(x_i-\mu)m_i^\mathsf{T}\}\{\sum_i E[Z_iZ_i^\mathsf{T}\mid x_i]\}^{-1}\). Equivalently, this is least squares with an intercept: subtract the sample averages of the observed vectors and the conditional factor means. The centered normal equations become \(\Lambda_{\mathrm{new}}A=B\), so \(\Lambda_{\mathrm{new}}=BA^{-1}\). The intercept equation is \(\mu_{\mathrm{new}}=\bar x-\Lambda_{\mathrm{new}}\bar m\). This establishes the joint updates rather than treating the mean as fixed during the loading maximization.

For the uniqueness update, let \(S_j=\sum_i E[(x_{ij}-\mu_j-\lambda_j^\mathsf{T}Z_i)^2\mid x_i]\). The part of the auxiliary function involving \(\psi_j\) is \(-\frac{n}{2}\log\psi_j-\frac{S_j}{2\psi_j}\). If \(S_j>0\), differentiating with respect to \(\psi_j>0\) and setting the derivative to zero gives the maximizer \(\psi_{j,\mathrm{new}}=S_j/n\). If \(S_j=0\), the auxiliary function increases without bound as \(\psi_j\) approaches zero, so there is no maximum over \(\psi_j>0\); a positive lower bound or another boundary treatment is needed. In computation, each expectation can be evaluated as \[ (x_{ij}-\mu_j-\lambda_j^\mathsf{T}m_i)^2+\lambda_j^\mathsf{T}V\lambda_j. \] The first term is the squared residual at the conditional factor mean; the second accounts for conditional factor uncertainty.

Why the Unknown-Mean Update Matters

The observed-data Gaussian likelihood separates into a mean component and a covariance component. This means that its maximum-likelihood estimate of the mean is \(\bar x\). But this fact does not make \(\bar x\) the general mean update in the EM M-step: the M-step maximizes the expected complete-data log-likelihood using the factor moments from the current E-step. The loading update and mean update must be compatible with those moments.

A convenient alternative is to center the data at \(\bar x\) before fitting the covariance model, then fix the mean at zero during EM. This uses the separation of the observed Gaussian likelihood into mean and covariance components. It is a different parameterization of the fitting calculation; it should not be confused with claiming that the unrestricted EM mean update equals \(\bar x\) at every iteration.

Likelihood Monotonicity and a Practical Pitfall

An exact EM iteration cannot decrease the observed-data likelihood, provided the E-step uses the current conditional distribution and the M-step does not decrease the auxiliary function. For this Gaussian model with positive definite \(\Psi\), the relevant densities and conditional expectations are well-defined. To see the likelihood argument, let \(q_i(z)\) be the old conditional density of \(Z_i\) given \(x_i\), and let \(p_\theta(x_i,z)\) be the joint density under parameter \(\theta\). The identity \[ \ell(\theta)-\ell(\theta_{\mathrm{old}}) =Q(\theta\mid\theta_{\mathrm{old}})-Q(\theta_{\mathrm{old}}\mid\theta_{\mathrm{old}}) +\sum_i \operatorname{KL}\!\left(q_i\,\middle\|\,p_\theta(\,\cdot\mid x_i)\right) \] follows by substituting the conditional-density identity \(p_\theta(z\mid x_i)=p_\theta(x_i,z)/p_\theta(x_i)\) into the definition of each Kullback–Leibler divergence. The divergence is nonnegative: for densities \(q,r\), Jensen's inequality gives \(E_q[\log(r/q)]\leq \log E_q[r/q]\leq 0\), so \(E_q[\log(q/r)]\geq0\). Thus an M-step that increases \(Q\) also does not decrease \(\ell\). Equality can occur, for example, at a stationary point. \(\square\)

• Common mistake: setting \(\mu_{\mathrm{new}}=\bar x\) automatically in the unrestricted EM M-step. That uses the observed-data mean estimate instead of maximizing the current auxiliary function jointly over \(\mu\) and \(\Lambda\). The correction is \(\mu_{\mathrm{new}}=\bar x-\Lambda_{\mathrm{new}}\bar m\), or to center the data and explicitly fix the mean at zero.
• Another common mistake: using \(m_im_i^\mathsf{T}\) as the conditional second moment. The correct quantity is \(V+m_im_i^\mathsf{T}\); omitting \(V\) discards uncertainty about the unobserved factors.
• Numerical caution: the updates are for a specified factor dimension and a diagonal, positive uniqueness matrix. Check that the resulting uniqueness estimates are admissible and monitor convergence; rotational nonuniqueness of the loading representation is not removed by these EM formulas.