Covariance Parameterization and Numerical Stability

Why Parameterize a Covariance Matrix?

In a generalized linear mixed model, random effects represent group-specific departures from the fixed-effects structure. Their covariance matrix determines how those departures vary and co-vary. Estimation is therefore not just a search for numerical values: every candidate covariance matrix must be mathematically admissible. In particular, a covariance matrix must be symmetric and positive semidefinite. A parameterization that produces invalid matrices can cause an optimizer to spend effort outside the model’s allowable parameter space.

The Cholesky parameterization addresses this constraint by expressing a positive-definite covariance matrix as a product of a lower-triangular matrix and its transpose. This gives the optimizer coordinates that generate valid matrices throughout the interior of the parameter space. It does not, by itself, guarantee that the likelihood is well behaved or that the estimated components are numerically stable. Those questions still require diagnostics.

For example, with a random intercept and a random slope, one group’s linear predictor may contain \(b_{0j}+b_{1j}x_{ij}\). The vector \(b_j=(b_{0j},b_{1j})^\mathsf{T}\) has a two-by-two covariance matrix. The response distribution is specified conditional on the random effects, and observations are ordinarily assumed conditionally independent given those effects and the covariates. The Cholesky factor concerns the distribution of the random effects; it does not change the chosen response distribution or link function.

How the Factor Guarantees Validity

For two random effects, write the lower-triangular factor as \(L\). Multiplying \(L\) by its transpose gives the covariance entries directly:

The matrix is symmetric by construction. Its first diagonal entry is \(a^2\), and its second is \(c^2+d^2\), so both are nonnegative. More importantly, its determinant is \(a^2d^2\). If \(a>0\) and \(d>0\), then both diagonal entries and the determinant are positive, which makes this two-by-two symmetric matrix positive definite. The same argument extends to \(q\) dimensions: for any nonzero vector \(v\),

because a triangular \(L\) with nonzero diagonal is invertible, so \(L^\mathsf{T}v\ne 0\). Thus the product representation guarantees positive definiteness when the factor has a nonzero diagonal. This is a validity result, not a guarantee that the data identify every covariance component well.

A frequent optimization convention is to set each diagonal entry to an exponential, such as \(a=\exp(\eta_1)\) and \(d=\exp(\eta_2)\), while leaving the off-diagonal entry \(c\) unconstrained. This allows the optimizer to work with unrestricted real-valued coordinates \(\eta_1,\eta_2,c\), while ensuring that \(a,d>0\). Implementations may use a scaled version of this construction or a different internal convention, so reported optimizer parameters need not equal the displayed entries of \(L\).

Cholesky Coordinates Are Not Correlations

The entries of \(L\) are computational coordinates. Except in special cases, they are not the standard deviations, variances, or correlations that researchers usually report. For a two-dimensional random effect with standard deviations \(\sigma_1,\sigma_2\) and correlation \(\rho\), one valid Cholesky factor is

Multiplication gives diagonal entries \(\sigma_1^2\) and \(\sigma_2^2\), and off-diagonal entry \(\rho\sigma_1\sigma_2\). Notice that \(L_{21}=\rho\sigma_2\), not \(\rho\), and \(L_{22}\) is not generally \(\sigma_2\). In dimensions above two, the relationship between a Cholesky factor and the resulting correlations is more involved still.

When fitting by maximum likelihood, re-expressing the covariance parameters in Cholesky coordinates does not add a Jacobian term to the likelihood. The likelihood is evaluated at the covariance matrix generated by the current optimizer coordinates. A Jacobian is relevant when transforming a probability density on parameters, as in a Bayesian prior transformation; that is a different calculation from maximizing the likelihood over a reparameterized model.

Scaling, Conditioning, and Diagnostics

A valid matrix can still be numerically awkward. Random-intercept and random-slope components often have different units and magnitudes, especially when the predictor is recorded on a large scale. The resulting covariance coordinates can differ greatly in size. This can make numerical optimization less well scaled and make small changes in one component difficult to compare with changes in another. Scaling the predictor changes the units and numerical values of the slope component; it does not change the underlying fitted linear predictor when the coefficients are transformed consistently.

Useful checks include confirming that the reported covariance matrix is symmetric, recalculating its standard deviations and correlations from its entries, and checking that its entries agree with the fitted factor under the software’s stated convention. Compare results across sensible predictor scales or equivalent parameterizations when numerical behavior is concerning. Also examine the convergence message, scaled gradient, Hessian, and optimizer sensitivity using the corresponding diagnostics. A convergence flag alone does not establish that every component is accurately estimated.

Choices and Limitations

Directly optimizing all covariance entries is easy to describe, but arbitrary entries need not produce a positive-semidefinite matrix. Constraints can be added, but they are themselves a numerical problem. A standard-deviation-and-correlation parameterization is often easier to interpret, but it also requires constraints: standard deviations must be nonnegative and the correlation structure must be valid. The Cholesky representation enforces positive definiteness through its factor structure and is consequently a common computational choice.

The choice is not cost-free. Cholesky entries depend on the ordering of the random effects, and their magnitudes depend on the units used. They are useful for computation but are usually less directly interpretable than variances, standard deviations, and correlations. Further, strictly positive diagonal entries represent positive-definite matrices in the interior; they do not represent every positive-semidefinite matrix exactly. A Cholesky parameterization therefore ensures validity of generated matrices, not adequate data support, good conditioning, or a scientifically appropriate random-effects structure.