Direct Likelihood Optimization and Numerical Failure Modes

Why Optimize the Observed-Data Likelihood Directly?

With incomplete multivariate data, each case contributes information about the coordinates that were observed, not a complete vector that can simply be analyzed as if nothing were missing. The Gaussian observed-data likelihood accounts for this by using the appropriate marginal normal distribution for each case’s observed coordinates. Directly optimizing that likelihood is an alternative to the EM algorithm: it can use general-purpose numerical optimizers and can be convenient when the model has additional constraints or parameters.

The numerical task is not just to find a parameter vector that makes an optimizer stop. Covariance matrices must remain positive definite, and the likelihood can be poorly behaved near the boundary of that parameter space. A small step in an unconstrained covariance entry can produce an invalid matrix; a nearly singular matrix can produce a very large likelihood or unreliable calculations. The objective, parameterization, and convergence checks therefore all matter.

Writing the Objective Correctly

For each case, let \(S_i\) be the coordinate-selection matrix that extracts the entries in \(O_i\). Write \(d_i=|O_i|\), \(\mu_i=S_i\mu\), and \(\Sigma_i=S_i\Sigma S_i^{T}\). The observed subvector has distribution \(N_{d_i}(\mu_i,\Sigma_i)\). Thus, up to the ignorable missingness mechanism, the log-likelihood is

The determinant and quadratic form are pattern-specific. In particular, a case with only one observed coordinate contributes the univariate marginal density for that coordinate; it does not contribute a complete-vector density with an invented value for the missing coordinate. Different cases can therefore involve different submatrices of \(\Sigma\). This is one reason a simple complete-data covariance formula generally does not maximize this objective.

Scores and Covariance-Safe Optimization

An optimizer may use analytic derivatives, numerical derivatives, or both. The following score expressions also provide a useful check on a numerical implementation. For each nonempty pattern, set \(r_i=Y_{i,O_i}-S_i\mu\) and \(A_i=(S_i\Sigma S_i^T)^{-1}\). Differentiating the log-likelihood gives

The second expression is the symmetric-matrix gradient under the trace inner product: its inner product with a symmetric perturbation \(H\) is the first-order change in log-likelihood in direction \(H\).

To verify the result, consider one nonempty pattern with \(C_i=S_i\Sigma S_i^T\), \(r_i=Y_{i,O_i}-S_i\mu\), and \(A_i=C_i^{-1}\). For a symmetric perturbation \(dC_i\), the determinant and inverse differentials are \(d\log|C_i|=\operatorname{tr}(A_i\,dC_i)\) and \(dA_i=-A_i(dC_i)A_i\). Consequently,

Rewriting the scalar quadratic differential as a trace shows that the gradient with respect to \(C_i\) is \(\tfrac12(A_i r_i r_i^T A_i-A_i)\). Since \(dC_i=S_i(d\Sigma)S_i^T\), the corresponding gradient with respect to \(\Sigma\) is \(S_i^T\{\tfrac12(A_i r_i r_i^T A_i-A_i)\}S_i\). For a mean perturbation \(d\mu\), \(dr_i=-S_i d\mu\), and the quadratic term has differential \(r_i^TA_iS_i d\mu\). Its gradient is \(S_i^TA_ir_i\). Summing these contributions over cases proves both formulas. \(\square\)

Optimizing the raw entries of \(\Sigma\) is hazardous: an unconstrained step can leave the positive-definite cone, making log-determinants and inverses undefined. A standard remedy is a Cholesky parameterization \(\Sigma=LL^T\), where \(L\) is lower triangular and its diagonal entries are strictly positive. For example, set \(L_{jj}=\exp(\eta_j)\), with unconstrained real parameters \(\eta_j\); the entries below the diagonal can also be optimized as unconstrained real numbers. Every finite parameter vector then gives a positive-definite covariance matrix.

Numerical Failure Modes and Diagnostics

A positive-definite parameterization prevents invalid covariance matrices, but it does not guarantee a well-behaved optimization problem. In particular, the likelihood can increase without bound as a covariance eigenvalue approaches zero when the data lie in a lower-dimensional configuration.

Poor initialization can also send an algorithm into a narrow, badly conditioned region. A covariance start with a very small eigenvalue makes inverses large and objective changes sensitive to tiny parameter movements. In a Cholesky parameterization, extremely negative diagonal log-parameters create the same practical difficulty even though the matrix remains positive definite. Multiple dispersed starts can reveal whether the optimizer is finding the same interior solution or repeatedly approaching different boundaries.

A small change in log-likelihood alone does not establish convergence to a useful maximum: a flat region, an ill-conditioned boundary, or a failed line search can also produce small changes. Likewise, a positive-definite estimate is not automatically well-conditioned. If a lower bound on eigenvalues, a penalty, or a covariance constraint is introduced to stabilize fitting, it changes the optimization problem and should be reported as part of the model or estimation procedure—not presented as though it were merely a computational detail.

Common Mistake

A common error is to treat an optimizer’s “converged” message as proof that a finite, reliable maximum-likelihood estimate has been found. That conclusion does not follow: an algorithm can stop near a singular boundary, at a point with a substantial gradient, or at a local solution whose quality depends on initialization. The correction is to examine the fitted covariance spectrum and condition number, check gradients and likelihood values, compare well-separated starts, and determine whether the likelihood appears to keep increasing as an eigenvalue shrinks. If it does, report the boundary or nonexistence issue rather than interpreting the returned finite iterate as an ordinary interior MLE.