Cholesky Factors for Positive-Definite Covariance Matrices

Why Factor a Random-Effects Covariance Matrix?

In an NLME model, a vector of subject-level random effects can represent how several individual parameters jointly depart from their typical values. The matrix \(\Omega\) describes their variances and covariances. When random effects are correlated, a valid covariance matrix must satisfy more than symmetry: it must also be positive definite if the model assumes a nonsingular multivariate normal distribution.

The Cholesky parameterization expresses \(\Omega\) as a product of a lower-triangular matrix and its transpose. This gives a direct way to build a valid covariance matrix from a set of factor entries. It also makes clear how correlated effects can arise from underlying independent standard normal variables.

For a vector of random effects \(\eta\), the usual covariance specification is \(\eta\sim N(0,\Omega)\), where the zero vector is the mean and \(\Omega\) is the covariance matrix. The factorization gives an equivalent construction: if \(z\sim N(0,I)\), then \(\eta=Lz\) has covariance \(\Omega\). Indeed, \(E(\eta)=0\) and

Thus, the factor does not remove correlation from the model. It constructs correlated random effects by combining independent standard normal components.

How the Factor Guarantees Positive Definiteness

Suppose \(L\) is lower triangular with positive diagonal entries. For any nonzero vector \(x\),

The last inequality holds because a triangular matrix with nonzero diagonal is invertible, so \(L^Tx\ne 0\) whenever \(x\ne 0\). Therefore \(LL^T\) is positive definite. This is the key computational advantage: choosing a lower-triangular factor with positive diagonal entries automatically produces a valid, nonsingular covariance matrix.

For a two-dimensional factor,

The diagonal entries of \(\Omega\) are variances, and its off-diagonal entries are covariances. The entries of \(L\), in contrast, are factor parameters; they are not generally variances, standard deviations, or correlations. For example, the first variance is \(l_{11}^2\), while the second includes contributions from both \(l_{21}\) and \(l_{22}\).

Recovering a Factor and Justifying the Construction

For a general \(p\)-dimensional positive-definite matrix, the entries of the factor can be found successively. For \(i\ge j\), the recursion is

where \(\omega_{ij}\) is the \((i,j)\) entry of \(\Omega\). The positive square root is used for every diagonal entry. Positive definiteness ensures each quantity under a square root is positive, so the recursion can proceed.

A brief block argument explains why the factorization exists for every positive-definite matrix. Assume the leading \((p-1)\)-dimensional block \(A\) has factorization \(A=RR^T\), and write the full matrix as \(\Omega=\begin{pmatrix}A&b\\b^T&d\end{pmatrix}\). Set \(v=R^{-1}b\). The scalar \(d-v^Tv=d-b^TA^{-1}b\) is positive: evaluate the positive quadratic form of \(\Omega\) at the nonzero vector \(\begin{pmatrix}-A^{-1}b\\1\end{pmatrix}\), which gives exactly \(d-b^TA^{-1}b>0\). Hence

has positive diagonal and direct multiplication gives \(LL^T=\Omega\). Starting with the one-dimensional case and applying this argument successively establishes existence. The recursion also determines each entry in order, and the positive diagonal requirement fixes each square-root choice, establishing uniqueness.

Using the Parameterization in NLME Models

An estimation procedure can work with the entries of \(L\) and compute \(\Omega=LL^T\) whenever it needs the covariance matrix. The off-diagonal entries of \(L\) may be positive, negative, or zero. The diagonal entries must remain positive; a common strategy is to represent each diagonal as the exponential of an unconstrained parameter. The exact implementation varies, but the purpose is the same: prevent the optimizer from proposing a factor that would produce an invalid covariance matrix.

If random effects have different units, the factor entries carry units implied by the products that form \(\Omega\). In common NLME specifications where the random effects are deviations on log scales, they are dimensionless. In either case, read the implied \(\Omega\) to interpret variances and covariances; do not interpret an off-diagonal entry of \(L\) as the correlation between two effects.

The factor depends on the ordering of the random effects. Reordering the effects generally changes \(L\), even when the reordered covariance matrix describes the same relationships. This is expected: \(L\) is a computational representation of \(\Omega\), not an order-free summary of the pairwise associations.

Diagnostics, Alternatives, and Limitations

A Cholesky parameterization ensures mathematical positive definiteness when its diagonal entries are positive, but it does not ensure that the data contain enough information to estimate every covariance reliably. Inspect the fitted covariance matrix, its standard deviations and correlations, and any convergence or numerical warnings. Very small diagonal elements, correlations close to \(1\) or \(-1\), large changes under small data or starting-value changes, and unstable estimates can indicate a nearly singular matrix or weak information about some combinations of effects.

A positive-definite matrix has strictly positive eigenvalues. If a computed matrix has a zero or negative eigenvalue, it is not strictly positive definite; rounding can also make a nearly singular matrix appear invalid. The Cholesky algorithm may fail or become numerically unstable when a pivot is nonpositive or very small. Such behavior should prompt investigation of numerical precision, the fitted model, and the information available in the data, rather than an assumption that the factorization itself has repaired the problem.

One alternative is to estimate variances and correlations directly. For two effects, this is straightforward if the correlation is constrained to lie between \(-1\) and \(1\). In higher dimensions, however, individually permissible pairwise correlations do not necessarily combine into a positive-definite correlation matrix. Directly estimating covariance entries likewise requires constraints to keep the complete matrix valid. The Cholesky approach enforces validity through the factor, though its entries may be less directly interpretable than standard deviations and correlations.