Covariance Parameterizations And Numerical Stability

Why the Covariance Parameterization Matters

In a shared random-effects joint model, the random effects connect an individual’s longitudinal trajectory to their event risk. Their covariance matrix describes both the variability of those effects and how they vary together. Estimation therefore requires not only finding plausible covariance values, but also keeping the matrix in the valid range throughout numerical optimization.

For example, a random intercept and random slope may be strongly related, but their covariance cannot be chosen independently of their variances. A proposed combination that violates the mathematical constraints does not describe any valid joint distribution. An optimizer that encounters such a combination may fail, return an unusable likelihood, or behave erratically. Choosing a parameterization that respects the constraints can prevent some failures, although it cannot resolve weak information in the data or guarantee successful estimation.

Definition: A covariance matrix \(D\) for a random-effects vector \(b_i\) is positive definite if \(v^\mathsf{T} D v > 0\) for every nonzero vector \(v\). It is positive semidefinite if \(v^\mathsf{T} D v \geq 0\) for every vector \(v\). A positive-semidefinite matrix may have zero eigenvalues; if it has any zero eigenvalues, it is singular and can represent a boundary case. Positive-definite matrices are also positive semidefinite.

Direct Covariance Estimation and Its Constraints

Let \(b_i=(b_{0i},b_{1i})^\mathsf{T}\) be a random intercept and random slope, with \(b_i\sim N(0,D)\). The direct covariance parameterization estimates the entries of \(D\) itself:

D = \begin{pmatrix} a & c \\ c & d \end{pmatrix}, \qquad a=\operatorname{Var}(b_{0i}),\quad d=\operatorname{Var}(b_{1i}),\quad c=\operatorname{Cov}(b_{0i},b_{1i}).

For this two-dimensional matrix to be positive definite, it must satisfy \(a>0\) and \(ad-c^2>0\). Equivalently, \(a>0\), \(d>0\), and the correlation \(c/\sqrt{ad}\) must be strictly between \(-1\) and \(1\). To allow a positive-semidefinite matrix, the conditions are \(a\geq0\), \(d\geq0\), and \(ad-c^2\geq0\). Directly optimizing \(a,c,d\) without constraints can step outside this feasible region.

In joint-model estimation, the random effects are typically integrated out of the likelihood. A candidate matrix that is not positive semidefinite cannot serve as a covariance matrix for the assumed random-effects distribution. Even when all iterates are valid, a matrix near singularity can make calculations sensitive: variation in one direction is nearly absent, and computations that rely on matrix inversion or a stable factorization can become poorly conditioned.

Worked Example: Checking a Directly Estimated Covariance Matrix

Suppose a fitted random-intercept/random-slope model reports \(D=\begin{pmatrix}0.64&0.12\\0.12&0.09\end{pmatrix}\), where the intercept is measured in outcome units and the slope in outcome units per year. Check whether this is a valid positive-definite covariance matrix and interpret its entries.

The diagonal entries are positive: \(a=0.64\) and \(d=0.09\). The determinant is \(ad-c^2=(0.64)(0.09)-(0.12)^2=0.0576-0.0144=0.0432\). Since this is greater than zero, the matrix is positive definite.

The random-intercept standard deviation is \(\sqrt{0.64}=0.80\) outcome units, and the random-slope standard deviation is \(\sqrt{0.09}=0.30\) outcome units per year. Their correlation is \(0.12/(0.80\times0.30)=0.50\), a positive association between subject-specific intercept and slope deviations. The covariance has units of outcome-squared per year; it is not itself a correlation.

Cholesky Parameterization

A Cholesky parameterization writes a positive-definite covariance matrix as \(D=LL^\mathsf{T}\), where \(L\) is lower triangular and has strictly positive diagonal entries. For two random effects, write \(L=\begin{pmatrix}\ell_{11}&0\\\ell_{21}&\ell_{22}\end{pmatrix}\). Then the covariance entries are determined by the factor:

D=LL^\mathsf{T} = \begin{pmatrix} \ell_{11}^2 & \ell_{11}\ell_{21}\\ \ell_{11}\ell_{21} & \ell_{21}^2+\ell_{22}^2 \end{pmatrix}.

To make the diagonal positive without imposing an explicit inequality constraint on each optimization step, an implementation can set \(\ell_{11}=\exp(\eta_1)\) and \(\ell_{22}=\exp(\eta_2)\), while \(\ell_{21}\) remains unconstrained. Every finite value of \((\eta_1,\eta_2,\ell_{21})\) then gives a positive-definite \(D\). This changes the optimization coordinates: it does not change which strictly positive-definite matrices can be represented.

Key result: If \(L\) is a square, nonsingular lower-triangular matrix, then \(D=LL^\mathsf{T}\) is positive definite. For any nonzero \(v\), \(v^\mathsf{T}Dv=v^\mathsf{T}LL^\mathsf{T}v=(L^\mathsf{T}v)^\mathsf{T}(L^\mathsf{T}v)>0\), because nonsingularity ensures \(L^\mathsf{T}v\ne0\).

This gives the central numerical advantage: an optimizer working over the Cholesky parameters cannot produce an indefinite covariance matrix through an ordinary finite update. The derivation also shows the limitation. If a diagonal factor approaches zero, the covariance approaches the boundary of the positive-semidefinite region. An exponential transform cannot attain a zero diagonal at any finite parameter value, so a boundary estimate may be approached only as a limit.

Worked Example: Converting Between a Covariance Matrix and Its Cholesky Factor

Use the positive-definite matrix from the previous example. Find a lower-triangular \(L\) such that \(D=LL^\mathsf{T}\), and verify the result.

The first diagonal element must satisfy \(\ell_{11}^2=0.64\), so choose \(\ell_{11}=0.80\). The off-diagonal entry requires \(\ell_{11}\ell_{21}=0.12\), giving \(\ell_{21}=0.12/0.80=0.15\). The lower-right entry then requires \(\ell_{21}^2+\ell_{22}^2=0.09\), so \(\ell_{22}^2=0.09-0.15^2=0.09-0.0225=0.0675\), and \(\ell_{22}=\sqrt{0.0675}\approx0.2598\).

Thus \(L=\begin{pmatrix}0.80&0\\0.15&0.2598\end{pmatrix}\), to four significant figures. Multiplying gives \(LL^\mathsf{T}\approx\begin{pmatrix}0.6400&0.1200\\0.1200&0.0900\end{pmatrix}\), with rounding accounting for any last-digit difference. In the exponential-diagonal form, \(\eta_1=\log(0.80)\approx-0.2231\) and \(\eta_2=\log(0.2598)\approx-1.348\).

How Invalid Updates Appear

A direct parameterization can be useful when the optimizer supports constrained optimization and handles the positive-definiteness conditions reliably. But if the entries are treated as unconstrained, a proposed update may be invalid even when the previous one was valid. The following calculation shows the issue:

Worked Example: Detecting an Indefinite Direct Update

Suppose an optimizer proposes \(D=\begin{pmatrix}0.04&0.025\\0.025&0.01\end{pmatrix}\). Determine whether this can be a covariance matrix.

Both variances are positive, but the determinant is \((0.04)(0.01)-(0.025)^2=0.0004-0.000625=-0.000225\). The determinant is negative, so one eigenvalue is negative and \(D\) is indefinite. The implied correlation would be \(0.025/(\sqrt{0.04}\sqrt{0.01})=0.025/(0.20\times0.10)=1.25\), which is outside the permissible range. This update cannot define a random-effects covariance matrix.

A Cholesky-based optimizer would instead update the factor entries and construct \(D=LL^\mathsf{T}\), preserving positive definiteness for finite positive diagonal factors. For instance, \(L=\begin{pmatrix}0.20&0\\0.05&0.08\end{pmatrix}\) gives \(D=\begin{pmatrix}0.0400&0.0100\\0.0100&0.0089\end{pmatrix}\). Its determinant is \(0.04(0.0089)-0.01^2=0.000256>0\), as expected.

Parameterization Choices, Diagnostics, and Limits

A common alternative is to parameterize standard deviations and correlations: \(D=SRS\), where \(S\) is diagonal with positive random-effects standard deviations and \(R\) is a correlation matrix. For two random effects this means estimating their standard deviations and a correlation restricted to \((-1,1)\). It is often easier to interpret than raw covariance entries, but it still requires a valid correlation parameterization. For more than two effects, independently restricting each pairwise correlation to \((-1,1)\) is not enough to guarantee that the full correlation matrix is positive definite. Cholesky parameterization directly addresses the full-matrix constraint.

The choice is not simply “Cholesky always works, direct estimation never does.” A constrained direct method can be stable and transparent; a Cholesky method guarantees validity of its constructed matrices but may still be poorly conditioned, difficult to optimize, or sensitive to parameter scaling. Different coordinates can also produce different optimization behavior even when they represent the same interior covariance matrices. A Cholesky factor is not usually interpreted as a set of variances and covariances; transform it back to \(D\) before reporting those quantities.

Inspect the reported covariance matrix, not only the optimizer’s convergence flag. Useful checks include:

Verify symmetry, positive diagonal entries, and positive eigenvalues for a claimed positive-definite estimate. For a two-dimensional matrix, check \(a>0\) and \(ad-c^2>0\).
Look for correlations extremely close to \(1\) or \(-1\), very small eigenvalues, or a large ratio of largest to smallest eigenvalue. These can indicate near-singularity, weakly determined directions, or scaling problems.
Check for non-finite entries, failed matrix factorizations, invalid likelihood evaluations, or abrupt changes in the objective when covariance parameters update.
Compare the implied variances, correlations, and units with the model’s random-effects definition. An impossible correlation or a variance in an unexpected scale can reveal an indexing or transformation error.

These symptoms do not diagnose one cause by themselves. Near-singularity may reflect a genuine boundary pattern in the data, weak identification, or numerical difficulty. A positive-definite parameterization prevents invalid interior iterates; it does not create information that the data do not contain, settle whether a boundary model is appropriate, or guarantee that all model components are identifiable. The earlier material on structural and practical identifiability is relevant when an apparently unstable covariance estimate may reflect limited information rather than an implementation defect.

Common pitfall: “The optimizer returned a covariance matrix, so it must be valid and reliably estimated.” This is not necessarily true. A reported matrix may be nearly singular, numerically rounded into an invalid form, or poorly identified. Check the matrix’s eigenvalues or positive-definiteness conditions and examine its scale and implied correlations; a valid matrix alone does not establish precise estimation.

Check Your Understanding

Answer each question, then compare with the solutions below.

For \(D=\begin{pmatrix}1&0.6\\0.6&0.25\end{pmatrix}\), is the matrix positive definite?
Explain why a finite Cholesky factor with strictly positive diagonal entries gives a positive-definite covariance matrix.
Can a positive-semidefinite matrix be nonsingular? When is it singular?
Why is restricting every pairwise correlation to \((-1,1)\) insufficient to ensure a valid correlation matrix when there are more than two random effects?
What does a very small eigenvalue of an estimated covariance matrix suggest, and what does it not prove?

Solutions. (1) The determinant is \(1(0.25)-0.6^2=0.25-0.36=-0.11\), so the matrix is not positive definite. (2) For nonzero \(v\), \(v^\mathsf{T}LL^\mathsf{T}v=\|L^\mathsf{T}v\|^2>0\), because a nonsingular \(L\) maps no nonzero vector to zero. (3) Yes. Positive-definite matrices are positive semidefinite and nonsingular; a positive-semidefinite matrix is singular when it has at least one zero eigenvalue. (4) Pairwise-valid correlations may combine into a matrix with a negative eigenvalue; the full matrix must satisfy positive-semidefiniteness, not merely pairwise bounds. (5) It suggests a nearly absent random-effects direction or numerical ill-conditioning. It does not by itself prove an implementation error or establish whether the cause is weak information, a boundary pattern, or poor scaling.