Parameterizing Correlation Matrices for Hamiltonian Monte Carlo

Why Correlations Need a Parameterization

In a GLMM with correlated random intercepts and slopes, the random effects for one group are often modeled as a multivariate normal vector. Their covariance matrix must be positive definite, and its correlations must be mutually compatible. HMC, however, updates variables in an unconstrained Euclidean space. Directly proposing each correlation independently does not guarantee a valid matrix: even if every proposed entry lies between \(-1\) and \(1\), the resulting matrix may fail to be positive definite.

The solution is to represent the correlation matrix through unconstrained coordinates and a smooth transformation. HMC then evolves those coordinates, while the transformation produces valid correlations at every finite point in the sampler's parameter space. The choice of coordinates also determines how a prior density must be expressed and which geometric difficulties may remain.

This tutorial concerns the parameterization of \(\Omega\), separately from the positive scale parameters in \(D\). The construction assumes a nonsingular correlation matrix in the interior of the valid space. Boundary cases, such as correlations exactly equal to \(1\) or \(-1\), are not attained at finite unconstrained coordinates. In Bayesian HMC, the log posterior in the unconstrained coordinates must be differentiable wherever the sampler explores it. As in other HMC applications, a transformation does not by itself guarantee good mixing.

A Two-Effect Transformation

For two random effects, there is only one correlation, \(\rho\), with \(-1<\rho<1\). A simple unconstrained coordinate is \(\eta\in\mathbb{R}\), transformed by the hyperbolic tangent:

Every finite \(\eta\) gives a valid correlation, and every interior correlation corresponds to \(\eta=\operatorname{atanh}(\rho)\). The derivative is positive, so the absolute Jacobian is

When the prior is specified as a density \(p_\rho(\rho)\) with respect to \(d\rho\), the density for \(\eta\) with respect to \(d\eta\) is \(p_\eta(\eta)=p_\rho(\tanh(\eta))\,|d\rho/d\eta|\). Equivalently, the log density used by HMC includes the log absolute Jacobian. This correction is required because a density changes when its coordinate scale changes.

Building a Higher-Dimensional Matrix from Partial Correlations

For more than two effects, a useful construction is based on sequential partial correlations. For \(q\) effects, define \(z_{ij}\in(-1,1)\), for \(i=2,\ldots,q\) and \(j=1,\ldots,i-1\). These coordinates are often called partial correlations in the Cholesky or hyperspherical construction. Transform each with \(z_{ij}=\tanh(\eta_{ij})\), where every \(\eta_{ij}\) is unconstrained.

Construct a lower-triangular matrix \(L\) row by row. For row \(i\), set

An empty product is 1, so \(L_{11}=1\). Each row has squared length 1: the terms are successive fractions of the remaining squared length, and the final diagonal term is the remainder. Consequently, \(\Omega=LL^\mathsf{T}\) has diagonal entries 1. Since \(L\) has positive diagonal entries for finite coordinates, it is nonsingular, so \(LL^\mathsf{T}\) is positive definite. Thus this construction produces a valid correlation matrix without constraints on the \(\eta_{ij}\).

Jacobian Corrections and Prior Densities

HMC evaluates a target density in the variables it updates. If the prior is defined on the correlation matrix, the density in unconstrained coordinates must include the Jacobian of the map. For the row-wise partial-correlation coordinates above, use the distinct off-diagonal entries of \(\Omega\) as coordinates on the correlation-matrix space. The absolute Jacobian from the \(z_{ij}\) coordinates to those entries is

This determinant can be obtained by ordering the correlation entries consistently with the sequential construction: each new partial correlation changes a new correlation by a remaining-length factor, giving a triangular derivative structure. For three effects, write \(a=z_{21}\), \(b=z_{31}\), and \(c=z_{32}\). Then \(\rho_{12}=a\), \(\rho_{13}=b\), and \(\rho_{23}=ab+\sqrt{1-a^2}\sqrt{1-b^2}c\). The determinant of the map \((a,b,c)\mapsto(\rho_{12},\rho_{13},\rho_{23})\) is \(\sqrt{1-a^2}\sqrt{1-b^2}\), matching the formula.

The further transformation \(z_{ij}=\tanh(\eta_{ij})\) contributes a factor \(1-z_{ij}^2\) for each coordinate. Therefore the complete Jacobian from \(\eta\) to the distinct correlations is

The log of this positive quantity is added to the log prior density when the prior is expressed with respect to correlation-matrix coordinates. If instead a prior is deliberately defined directly on the \(\eta_{ij}\), its density is already on the HMC scale and this change-of-variables correction is not added. Those are different prior specifications, not interchangeable computational choices.

Prior Choice, Diagnostics, and Alternatives

A common prior for a correlation matrix is the LKJ family, whose density with respect to the distinct correlations is proportional to \(|\Omega|^{\xi-1}\), where \(\xi>0\). The determinant is a measure of how far the matrix is from singularity. When \(\xi=1\), the density is uniform over the valid correlation-matrix space; it does not mean that each pairwise correlation is independently uniform. In unconstrained coordinates, the LKJ log density and the transformation log Jacobian must both be included if the prior was specified on \(\Omega\).

For this construction, \(|\Omega|=\prod_{i>j}(1-z_{ij}^2)\). Combining this determinant with the Jacobian shows how an LKJ prior induces marginal structure on the partial-correlation coordinates. In particular, a simple prior on each \(\eta_{ij}\) is generally not equivalent to an LKJ prior. Prior interpretation therefore requires stating the scale on which the prior is defined.

Practical checks should include verifying that the transformed matrix is symmetric with unit diagonal, that its Cholesky diagonal is positive, and that the log Jacobian is included exactly once when required. Extreme values of \(|\eta_{ij}|\) make \(\tanh(\eta_{ij})\) very close to a boundary, where derivatives become small and floating-point calculations can lose precision. Evaluate expressions such as \(\log(1-z_{ij}^2)\) with numerically stable routines when available, and inspect the resulting matrices and posterior draws for near-singularity.

A frequent mistake is to say, “Each correlation is between \(-1\) and \(1\), so the full matrix is valid.” That condition is necessary but not sufficient: a collection of individually admissible correlations can yield a non-positive-definite matrix. The correction is to construct the matrix through a parameterization that guarantees positive definiteness, or explicitly constrain the full matrix as a valid correlation matrix.

Directly sampling constrained correlation entries can be conceptually simple, but boundaries and positive-definiteness constraints are awkward for unconstrained HMC. A Cholesky-based construction gives smooth unconstrained coordinates and a straightforward validity guarantee, but can have geometry that is difficult for HMC, especially near singular matrices. Other valid constructions are possible; whichever is chosen, the model's prior must be transformed consistently. No parameterization removes the need to assess posterior computation, and transformation alone does not resolve weak information about correlations or dependence between correlations and the random-effect scales.