Maintaining Positive-Definite Random-Effect Covariance

Why Covariance Validity Matters

In a nonlinear mixed-effects model, the random-effect covariance describes both the variability of individual effects and their dependence. For a \(q\)-dimensional random effect \(b_i\), the usual Gaussian specification is \(b_i\sim N_q(0,\Omega)\), where \(\Omega\) must be symmetric and positive definite. An update that violates either property is not merely an awkward numerical result: it is not a valid covariance matrix for this model. It can also make likelihood calculations fail, since these typically require \(\Omega^{-1}\) and \(\log|\Omega|\).

The issue in SAEM is subtle. A matrix of accumulated second moments can be positive semidefinite without being positive definite. Furthermore, averaging valid matrices preserves positive semidefiniteness only when the starting matrix is also positive semidefinite and the update uses a common convex weight. Separate updates to matrix entries do not have that guarantee. A reliable implementation must distinguish these facts and ensure that the covariance used in the model remains admissible.

Model and SAEM Covariance Statistic

Assume each subject has a \(q\)-dimensional random effect \(b_i\sim N_q(0,\Omega)\), with \(\Omega\) symmetric positive definite. In the covariance part of the complete-data log likelihood, terms involving \(\Omega\) are proportional to

Here \(n\) is the number of subjects. SAEM replaces the unavailable complete-data second moments by a stochastic approximation of the simulated sufficient statistic. Write \(A_k\) for the normalized statistic supplied at iteration \(k\), such as \(n^{-1}\sum_i b_i^{(k)}b_i^{(k)\mathsf T}\) from one simulated draw per subject, and \(S_k\) for its running approximation. The matrix covariance update is then

Each outer product \(b_i b_i^{\mathsf T}\) is positive semidefinite because \(v^{\mathsf T}b_i b_i^{\mathsf T}v=(b_i^{\mathsf T}v)^2\geq 0\). Thus averages of these outer products are positive semidefinite. They need not be positive definite: for example, a single draw in two dimensions gives a matrix of rank at most one. When \(S_k\) is positive definite, the covariance M-step sets \(\Omega_k=S_k\).

Why a Matrix Update Can Fail

A covariance matrix is constrained as a whole. Its entries cannot be chosen independently: for a two-dimensional covariance, the squared covariance must not exceed the product of the two variances. Updating diagonal entries and an off-diagonal entry with unrelated gains can break this constraint even if both the old matrix and the new matrix of proposed entries were intended to represent covariances.

When Matrix Averaging Preserves Validity

The relevant guarantee is about positive semidefiniteness, and it has an essential initialization condition. If \(S_{k-1}\) and \(A_k\) are both positive semidefinite and \(0\leq\gamma_k\leq1\), then their common-weight update remains positive semidefinite. It is not enough that \(A_k\) is positive semidefinite: the starting statistic must be as well.

For the first claim, take any vector \(v\). By the assumptions,

The result follows by induction from \(S_0\). For the stronger claim, if \(v\neq0\), then \(v^{\mathsf T}S_{k-1}v>0\). Since \(1-\gamma_k>0\) and the other quadratic form is nonnegative, \(v^{\mathsf T}S_kv>0\). At \(\gamma_k=1\), the update equals \(A_k\), which may be singular. This distinction explains why a positive-semidefinite SA statistic does not automatically give a positive-definite covariance at every iteration.

Cholesky Parameterization

A direct way to represent a positive-definite covariance is through its Cholesky factor. Every symmetric positive-definite \(\Omega\) has a unique factorization \(\Omega=LL^{\mathsf T}\), where \(L\) is lower triangular and has strictly positive diagonal entries. Conversely, every such \(L\) gives a positive-definite matrix: for nonzero \(v\), \(v^{\mathsf T}LL^{\mathsf T}v=\|L^{\mathsf T}v\|^2>0\), since \(L\) is invertible.

For unconstrained optimization over factor entries, write each diagonal entry as an exponential, \(L_{jj}=\exp(\eta_j)\), and leave the entries below the diagonal unconstrained. Any finite \(\eta_j\) gives a strictly positive diagonal entry. The transformed objective can then be optimized in the unconstrained coordinates \(\eta_j\) and \(L_{ij}\) for \(i>j\), while \(\Omega=LL^{\mathsf T}\) stays positive definite at each evaluation.

Using the Parameterization in SAEM

There are two related but distinct operations. SAEM commonly smooths sufficient statistics using the matrix recursion above, then performs an M-step. When the smoothed second-moment statistic is positive definite, the Gaussian covariance M-step is its value. Alternatively, an implementation can optimize the covariance-dependent objective in Cholesky coordinates, ensuring that every trial covariance is positive definite. Reparameterizing the M-step does not mean that one can average Cholesky factors with the usual matrix-statistic gain and obtain the same SAEM update: the map \(L\mapsto LL^{\mathsf T}\) is nonlinear, and factor averaging does not generally equal covariance averaging.

A short justification for the covariance M-step clarifies the boundary case. Let \(S=n^{-1}\sum_i b_i b_i^{\mathsf T}\) be the normalized second moment. Apart from a positive multiplier, maximizing the Gaussian objective is equivalent to maximizing \(-\log|\Omega|-\operatorname{tr}(\Omega^{-1}S)\). Its matrix derivative with respect to symmetric \(\Omega\) is \(-\Omega^{-1}+\Omega^{-1}S\Omega^{-1}\). Setting this to zero and multiplying on the left and right by \(\Omega\) gives \(S=\Omega\). This is an interior positive-definite solution when \(S\) is positive definite. If \(S\) is singular, that interior solution is unavailable; forcing a positive-definite answer requires a constrained or regularized procedure, not simply declaring the singular moment matrix positive definite.

One common safeguard is a lower bound on Cholesky diagonal entries, or an eigenvalue floor such as replacing an estimated matrix's eigenvalues below a chosen threshold with that threshold. Another is to use \(\Omega=LL^{\mathsf T}+\epsilon I\) for fixed \(\epsilon>0\). Each enforces positive definiteness, but imposes a constraint or changes the covariance estimate; it is not an exact unconstrained Gaussian M-step. The threshold should be chosen in relation to parameter scales, reported, and checked for sensitivity.

Diagnostics, Alternatives, and Limitations

Check covariance validity at the point where it is constructed and before computing its inverse or log determinant. Useful diagnostics include symmetry error, the smallest eigenvalue, Cholesky-factorization success, and the condition number. A matrix can be positive definite but nearly singular, creating unstable inverses and unreliable numerical optimization. Conversely, a small negative eigenvalue on the order of floating-point error may reflect rounding; symmetrizing with \((\Omega+\Omega^{\mathsf T})/2\) addresses asymmetry but does not by itself ensure positive definiteness.

If an update fails, trace the complete operation rather than clipping a reported correlation alone. Confirm that the sufficient statistic was formed from outer products consistently, the initial statistic was valid, a single gain was applied to the whole matrix, and no entrywise rounding or scaling changed the matrix. When a Cholesky factorization fails, inspect the eigenvalues and units; an arbitrary diagonal addition may conceal a coding error or materially change a variance estimate. Cholesky parameterization also cannot eliminate weak information in the data: it ensures admissibility, not accurate estimation, good mixing, or well-conditioned curvature.

Common Mistake