Structural Zeros in Correlation Matrices: A Thirteen-Parameterization Bake-Off

1. The Problem with Zeros in Correlation Matrices

Sometimes you know, before seeing any data, that certain correlations are exactly zero. For example, a priori knowledge that two variables are marginally¹ unrelated occurs in graphical models as either an explicit assumption or by pre-existing knowledge that rules out the variable coupling. Building this structure directly through the correlation matrix parameterization is what this blog post is about.

Sampling from an LKJ prior and setting selected cells to zero does not work because the resulting matrix will generally not be positive definite. Instead, the zeros have to be part of the parameterization. In this post I compare thirteen Stan (Carpenter et al. 2017) implementations that target the same distribution. I check that they agree numerically, then compare divergences, effective sample size, gradient evaluations, and run time.

The implementations fall into two groups. The first imposes the zero equalities sequentially in Cholesky or C-vine coordinates. The second samples each row in the subspace satisfying all of its zero constraints. I also introduce a triangular basis for the second approach. In these examples it removes the divergences seen in the entry-wise construction without the loss of efficiency I found with a QR basis.

2. Family 1: Forced Zeros

Seven members of the first group use the constrained-Cholesky approach I introduced in Pinkney (2024a) and presented at StanCon (Pinkney 2024b). The original method also allows a correlation to be fixed at a nonzero value or restricted to an interval. Writing \(\Omega = L L'\) with unit-norm rows of \(L\), the correlation in cell \((i,j)\) is \[ \Omega_{ij} = b_1 + L_{ij} L_{jj}, \qquad b_1 = \sum_{k<j} L_{ik} L_{jk}, \] This is affine in \(L_{ij}\), so setting the correlation to zero gives \(L_{ij}=-b_1/L_{jj}\) and uses no parameter. A bound on \(\Omega_{ij}\) similarly becomes a bound on \(L_{ij}\). Conditioning on \(\Omega_{ij}=0\) contributes \(-\log L_{jj}\) to the log density for each zero. The carpenter variant below is Bob Carpenter’s compact formulation of this case from a Stan Discourse discussion:

if (is_zero(i, j, zeros, zero_idx)) {
  real b = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
  L[i, j] = -b / L[j, j];   // implies Omega[i, j] == 0
  jacobian += -log(L[j, j]);
  zero_idx += 1;
}

The free entries can be parameterized in several ways. I compare the following nine variants:

label	free-cell parameterization
carpenter	direct Cholesky entries, bounded by the remaining stick
pinkney	correlation-scale value divided by \(L_{jj}\), multiplicative stick (Pinkney 2024a, 2024b)
ldl / ldl2 / ldl3	\(LDL'\) factorization (Pinkney 2024a); stick tracked via log_diff_exp, a log1m ratio, or the natural scale
cvine	C-vine partial correlations with the LKJ Beta factorization (Joe 2006; Lewandowski et al. 2009)
cvine_lookahead	exact equality coordinates, all-future 3-by-3 look-ahead, and a reduced free-cell Jacobian
rowscale	direct Cholesky entries, row-scaling stick
schur	canonical partial correlations, \(L_{ij} = r \tanh(x)\)

The seven Cholesky and \(LDL'\) variants share the same constraint. The value \(-b_1/L_{jj}\) is computed without reference to the remaining row norm, \(\sqrt{1 - \sum_{k<j} L_{ik}^2}\), and can exceed the available range. A later square root then receives a negative argument. In the \(K=7\) benchmarks below, these variants have divergences in about 1% of post-warmup transitions.

The cvine_lookahead model instead represents every zero as an exact ordinary-correlation equality, omits those coordinates from the parameter vector, and uses exact \(3\times3\) elliptope projections to narrow each free partial-correlation interval before sampling it. Its Jacobian contains only the free ordinary-correlation coordinates, and the LKJ kernel is added through \((\eta-1)\log\det(\Omega)\). The look-ahead substantially reduces infeasible proposals, although it does not remove every nonlocal feasibility boundary for dense patterns.

The three LDL variants are also useful for judging Monte Carlo variability. They define the same density and differ only in how some floating-point operations are arranged, but their total divergences over five patterns range from 356 to 466. Differences of that size should not be interpreted as evidence that one parameterization is better from a single set of fits.

3. Family 2: Exact-Constraint Subspaces

The second approach comes from Seth Axen’s Turing.jl discussion of positive definite matrices with structural zeros (Axen 2023):

For each column j, each column i for i < j that is orthogonal to column j decreases the dimensionality of the hemisphere by 1. So we simply construct a point on a lower dimensional hemisphere and use the QR decomposition of the orthogonal columns to lift it to a higher dimensional hemisphere.

Row \(i\) of a correlation Cholesky factor is a unit vector with a positive diagonal element. Since row \(j\) is zero after column \(j\), the condition \(\Omega_{ij}=0\) is equivalent to requiring rows \(i\) and \(j\) to be orthogonal. If row \(i\) has \(n_z\) constrained correlations, it can be sampled on the unit hemisphere in the orthogonal complement of those \(n_z\) rows. The dimension of this hemisphere is \(d=i-1-n_z\). The joint conditioning term for the row is \(-\tfrac12\log\det(A'A)\), where the columns of \(A\) are the rows to which it must be orthogonal.

My first two implementations, qr and qr2, construct this complement with a Householder QR decomposition. They use either a gnomonic map or the hyper-tanh peel (Pinkney 2026) for the hemisphere coordinates. Both models have zero divergences in the \(K=7\) benchmark, but their minimum ESS is about 40% of the entry-wise models. Changing the hemisphere map has little effect. The larger difference comes from the complement basis: the trailing columns of \(Q\) change as the preceding rows change, so the interpretation of the unconstrained coordinates also changes during sampling.

4. A Canonical Triangular Basis

I instead construct a basis tied to the free matrix entries. For each free column \(k\) of row \(i\), I start with the coordinate vector \(e_k\) and solve for the constrained columns by back-substitution. This is the entry-wise construction from the previous section applied to a basis vector rather than to the sampled row. If \(Z\) and \(F\) denote the constrained and free columns, the resulting basis matrix \(T\) satisfies \(T_F=I\) and \(T_Z=-(A_Z')^{-1}A_F'\).

To orthonormalize this basis, let \(S=\operatorname{chol}(T'T)\). Then \(v=T(S')^{-1}z\) preserves the norm of \(z\). I parameterize \(z\) using the row transformation from Pinkney (2025). With \(t=y'y\), \[ z = \frac{\sqrt{t+2}}{t+1}y, \qquad L_{ii}=\frac{1}{t+1}, \] and the log Jacobian is \(\tfrac12(d-2)\log(t+2)-(d+1)\log(1+t)\).

Two identities reduce the amount of linear algebra required.

1. The zero-partner rows have triangular support, so \(A_Z\) is triangular with the partners’ diagonals on its diagonal. Sylvester’s determinant identity gives \[ \det(A'A) = \Big(\prod_{j \in Z} L_{jj}\Big)^{2} \det(T'T). \] Thus \(-\tfrac12\log\det(A'A)\) becomes \(-\sum_j\log L_{jj}-\sum_k\log S_{kk}\). The second term is available from the Cholesky factor already used to orthonormalize \(T\).

2. Since \(T_F=I\), \(T'T=I+W'W\), where \(W\) contains only the \(n_z\) constrained rows. The cross-product therefore uses vectors of length \(n_z\) rather than \(i-1\).

Each constrained row requires one \(d\times d\) cholesky_decompose and several matrix-vector products. I call this model tri; the Stan code is in the appendix.

I also test an over-parameterized version, called donut, that uses one additional parameter per row. It normalizes a free vector, \(u=y/\lVert y\rVert\), before mapping it into the constrained subspace. The normalization is rotationally symmetric but singular at the origin. Following Axen (2022), I replace the usual radial term \(-\tfrac12\lVert y\rVert^2\) with \(-\tfrac12(\lVert y\rVert-m)^2\) so that the typical radius is away from zero. Since this term only changes the radius, it does not change the induced distribution on \(\Omega\). The largest difference in a free-cell posterior mean relative to carpenter is 0.0046. Across the five \(K=7\), \(\eta=4\) patterns:

radial prior	total divergences	median min ESS	median ESS / 1k gradients
\(m = 0\) (implicit chi)	2,549	4,778	19
\(m = 2\)	349	8,106	80
\(m = 4\)	8	10,778	192
\(m = 6\)	0	9,045	162

Increasing \(m\) removes the divergences and improves efficiency up to \(m=4\); \(m=6\) gives zero divergences in these runs. I use \(m=6\) below. The tri model does not require this additional choice and uses \(K-1\) fewer parameters.

5. The Bake-Off at K = 7

Setup: \(K = 7\), \(\eta = 4\), prior-only, five random zero patterns with roughly half of the 21 lower-triangular cells zeroed, 4 chains \(\times\) 2000 draws after 1000 warmup, default adapt_delta = 0.8. Correctness was checked on a held-out pattern first: all thirteen posteriors match the carpenter baseline within Monte Carlo error, and the constrained cells come back as numerically exact zeros (\(\sim 10^{-19}\)).

Divergences out of 8,000 post-warmup transitions:

pattern	carpenter	pinkney	ldl	ldl2	cvine	cvine_lookahead	rowscale	ldl3	schur	qr	qr2	tri	donut
seed 1	34	53	53	26	41	0	44	40	40	0	0	0	0
seed 2	61	65	74	107	64	0	78	104	60	0	0	0	0
seed 3	77	44	67	57	69	1	44	68	59	0	0	0	0
seed 4	137	157	130	112	186	0	155	152	123	0	0	0	0
seed 5	110	94	114	54	141	7	114	102	29	0	0	0	0
total	419	413	438	356	501	8	435	466	311	0	0	0	0

Min bulk ESS over the free cells of \(\Omega\) (8,000 draws):

pattern	carpenter	pinkney	ldl	ldl2	cvine	cvine_lookahead	rowscale	ldl3	schur	qr	qr2	tri	donut
seed 1	11639	10566	9979	10618	10219	11572	11319	9840	11036	3995	3362	11238	10676
seed 2	10817	9866	10826	9772	10188	12308	10153	10361	10040	4049	3728	11923	9045
seed 3	9344	11014	10244	10314	10946	9535	9420	10625	10787	4334	4304	10655	9429
seed 4	8395	9159	8326	9257	8771	9224	9463	8266	9458	3701	3815	9400	8161
seed 5	10381	10010	9268	8760	10301	10454	9218	9896	4770	4267	4153	11845	8619
median	10381	10010	9979	9772	10219	10454	9463	9896	10040	4049	3815	11238	9045

Efficiency (medians over the five patterns):

	carpenter	pinkney	ldl	ldl2	cvine	cvine_lookahead	rowscale	ldl3	schur	qr	qr2	tri	donut
min ESS / 1000 gradients	188.5	179.6	185.9	182.2	188.9	194.1	176.3	182.1	182.0	75.0	72.5	212.3	161.8
total wall clock, 5 fits (s)	4.53	4.10	4.40	4.51	5.14	2.65	4.21	4.51	4.15	2.07	2.06	1.86	2.39

The cvine_lookahead model records 8 divergences across the five patterns, compared with 501 for cvine, and completes the five fits in 2.65 rather than 5.14 seconds. The tri model still has zero divergences for all five patterns, the highest median minimum ESS, the highest median ESS per 1,000 gradients, and the shortest total run time. The QR models also remove divergences, but their changing complement basis produces substantially lower ESS. The triangular basis avoids that loss in this benchmark.

6. What the Look-Ahead Changes: cvine vs cvine_lookahead

The original cvine samples the unconstrained C-vine partial correlations from their LKJ Beta factors. When an ordinary correlation is fixed to zero, it solves the affine partial-correlation recursion for the required value. Nothing restricts earlier free partials to ensure that the solution will remain in \((-1,1)\), so the transform can meet a feasibility boundary late in the recursion.

cvine_lookahead carries intervals for the unsampled conditional correlations. Before selecting the current partial correlation \(x\), it uses the exact \(3\times3\) condition

\[ 1 + 2xyr - x^2 - y^2 - r^2 \ge 0 \]

to project every later pair onto a feasible interval for \(x\). After selecting \(x\), it tightens the later intervals and continues. Fixed ordinary correlations are point intervals, so zeros take no raw parameter in either model. The new transform uses the reduced Jacobian from raw values to the free ordinary correlations and adds the LKJ kernel as \((\eta-1)\log\det(\Omega)\); it therefore targets the same conditional LKJ distribution as cvine.

\(K=7\), \(\eta=4\), five patterns	cvine	cvine_lookahead
total divergences	501	8
median minimum bulk ESS	10,219	10,454
median ESS / 1,000 gradients	188.9	194.1
total wall clock (s)	5.14	2.65

The scaling sweep tells the same story. The original initializes in 12 of 20 configurations and fails for every \(K=15\) and \(K=25\) fit. cvine_lookahead initializes in 16 of 20 and reaches \(K=15\). On the 12 configurations where both initialize, divergences fall from 20,842 to 993, median minimum ESS rises from 6,447 to 7,876, and median ESS per 1,000 gradients rises from 84.3 to 152.4. It still fails under default random initialization at \(K=25\): enforcing every local \(3\times3\) condition is necessary, but for \(K\ge4\) it is not sufficient for a globally feasible completion satisfying all remaining zero constraints.

7. How Far Can \(3\times3\) Look-Ahead Propagate?

The gist implementation already has maximum breadth: at each C-vine pair \((i,j)\) it checks every future \(k>j\), not only \(k=j+1\). What it does not do is repeatedly propagate every narrowing through other triangles. A bound narrowed in one triangle can narrow another edge in a second triangle, which can then feed a third. To test that depth, I added an experimental transform that sweeps all currently active triples symmetrically for a fixed number of passes before each sequential coordinate. Each synchronous pass reads the previous pass’s intervals, intersects all three projected edges, and never moves an equality or an already sampled coordinate.

The bounds still depend only on earlier raw coordinates and unsampled boxes. The transformation remains triangular, so the Jacobian needs only the current interval width, logistic derivative, and the existing C-vine scale factor; there is no additional “projection Jacobian.” Fixed pass counts also avoid a tolerance-dependent stopping rule inside HMC.

The next table uses the same five \(K=7\), \(\eta=4\) patterns and sampler settings. Spans 1 through 5 limit both the pre-selection intersection and post-selection tightening to that many future indices. all reaches the maximum possible distance of 6 and is cvine_lookahead. The sweep rows retain all-future breadth and add synchronous all-triple propagation before every coordinate. For rowscale, pass 0 and pass 1 are the same executable with propagation disabled or enabled, which controls for floating-point changes between Stan programs.

experiment	total divergences	median min ESS	median ESS / 1k gradients	total time (s)
C-vine, span 1	206	11,153	214.5	3.80
C-vine, span 2	43	10,562	205.7	3.06
C-vine, span 3	9	10,454	194.1	2.93
C-vine, span 4	8	10,454	194.1	2.84
C-vine, span 5	8	10,454	194.1	3.05
C-vine, all future \(k\)	8	10,454	194.1	2.99
C-vine, all + 1 sweep	7	11,255	219.4	19.89
C-vine, all + 2 sweeps	7	11,255	219.4	37.96
C-vine, all + 4 sweeps	7	11,255	219.4	74.42
rowscale, 0 sweeps	362	9,567	174.8	4.54
rowscale, 1 sweep	207	10,334	188.6	9.11
rowscale, 2 sweeps	207	10,334	188.6	16.67
rowscale, 4 sweeps	207	10,334	188.6	31.39

Breadth is the useful part for the C-vine: nearest-neighbor look-ahead leaves 206 divergences, two-step leaves 43, three-step leaves 9, and spans 4, 5, and all future indices agree at 8. One synchronous all-triple pass before each coordinate changes the C-vine path slightly and removes one more divergence, but costs nearly seven times as much; passes 2 and 4 produce exactly the same diagnostics and ESS as pass 1 while only multiplying run time. The naive repeated sweep is \(O(pK^5)\) for \(p\) passes, so it is an experiment rather than a recommended default.

The projection also transfers to an entry-wise Cholesky transform, provided it is applied to ordinary correlations rather than to Cholesky entries. For rowscale, I maintain boxes for \(\Omega_{ij}\), sweep all same-row pairs before each coordinate against the already known \(\Omega_{jk}\), then map the contracted interval back through \(\Omega_{ij}=b_1+L_{ij}L_{jj}\). One pass cuts divergences from 362 to 207 and improves median efficiency; additional passes saturate immediately. It helps, but unlike the subspace construction it does not remove the nonlocal constrained-completion boundary.

I then ported the same one-pass ordinary-correlation sweep to six entry-wise transforms from the bake-off. At a common parity-check point, pass-zero mode reproduced each original log density and full raw-gradient vector to printed precision. The table compares pass zero with pass one in those same executables; this is a cleaner intervention than comparing separately compiled programs. The contraction happens on \(\Omega\), after which each transform retains its existing Cholesky or \(LDL'\) scale and equality Jacobians.

model	successful patterns, 0 \(\to\) 1 pass	divergences, 0 \(\to\) 1 pass	median ESS / 1k gradients, 0 \(\to\) 1 pass
carpenter	5 \(\to\) 5	398 \(\to\) 195	199.3 \(\to\) 202.6
pinkney	5 \(\to\) 5	395 \(\to\) 322	192.3 \(\to\) 199.1
ldl2	5 \(\to\) 5	389 \(\to\) 275	187.7 \(\to\) 190.3
ldl3	5 \(\to\) 5	370 \(\to\) 230	182.1 \(\to\) 190.9
rowscale	5 \(\to\) 5	362 \(\to\) 207	174.8 \(\to\) 188.6
schur	4 \(\to\) 5	224 \(\to\) 1932	191.9 \(\to\) 202.2

Across the 29 matched model-pattern pairs, one sweep reduces divergences from 2,138 to 1,422, raises median minimum ESS from 10,222 to 10,334, and raises median ESS per 1,000 gradients from 189.7 to 190.9. Wall time doubles from 25.9 to 52.1 seconds. The direction is favorable for every model in aggregate, but the gain is much smaller than for cvine_lookahead, and individual fits can move the other way because HMC trajectories are chaotic near the boundary. Local \(3\times3\) pruning is therefore useful as a guardrail for the entry-wise family, not a replacement for a parameterization that is feasible by construction.

On a held-out \(K=5\) pattern, the six one-pass posteriors also agree with the carpenter baseline: the largest free-cell mean difference is 0.016, the largest difference is 2.85 combined Monte Carlo standard errors, and the constrained cells remain zero to about \(10^{-19}\). The experimental Stan sources, benchmark scripts, and raw result CSVs are stored alongside this post.

8. Does It Hold Up? Scaling K and \(\eta\)

I next repeat the comparison for \(K\in\{5,7,9,15,25\}\) and \(\eta\in\{1,2,3,4\}\). This sweep includes five entry-wise models, both C-vine models, qr2, tri, and the donut model with \(m=6\). For each \(K\), I generate one Bernoulli(0.5) zero pattern and reuse it for all four values of \(\eta\). The realized number of zeros ranges from 6 of 10 cells at \(K=5\) to 143 of 300 at \(K=25\). Sampler settings are unchanged. There is one fit per configuration, so small differences should be interpreted with the variability from Section 2 in mind.

library(ggplot2)
library(dplyr)
library(tidyr)

source("../../R/theme_blog.R")

sweep <- read.csv("sweep-results.csv") |>
  mutate(
    family = case_when(
      model == "tri" ~ "Subspace (triangular basis)",
      model == "qr2" ~ "Subspace (QR basis)",
      model == "donut" ~ "Subspace (donut)",
      model == "cvine_lookahead" ~ "C-vine (look-ahead)",
      .default = "Forced zeros"
    ),
    div_pct = 100 * divergences / 8000,
    eta_lab = factor(paste0("eta == ", eta),
                     levels = paste0("eta == ", sort(unique(eta))))
  )

fam_colors <- c(
  "Forced zeros" = blog_colors$gold,
  "C-vine (look-ahead)" = blog_colors$red,
  "Subspace (QR basis)" = blog_colors$purple,
  "Subspace (donut)" = blog_colors$blue,
  "Subspace (triangular basis)" = blog_colors$teal
)

ggplot(sweep |> filter(!failed),
       aes(x = K, y = div_pct, group = model, color = family)) +
  geom_line(aes(linewidth = family, alpha = family, linetype = family)) +
  geom_point(size = 1.6) +
  facet_wrap(~eta_lab, nrow = 1, labeller = label_parsed) +
  scale_color_manual(values = fam_colors, name = NULL) +
  scale_linewidth_manual(guide = "none", values = c(
    "Forced zeros" = 0.6, "C-vine (look-ahead)" = 0.9, "Subspace (QR basis)" = 1.1,
    "Subspace (donut)" = 1.1, "Subspace (triangular basis)" = 1.3)) +
  scale_alpha_manual(guide = "none", values = c(
    "Forced zeros" = 0.65, "C-vine (look-ahead)" = 1, "Subspace (QR basis)" = 1,
    "Subspace (donut)" = 1, "Subspace (triangular basis)" = 1)) +
  scale_linetype_manual(guide = "none", values = c(
    "Forced zeros" = "solid", "C-vine (look-ahead)" = "solid", "Subspace (QR basis)" = "solid",
    "Subspace (donut)" = "22", "Subspace (triangular basis)" = "solid")) +
  scale_x_continuous(breaks = c(5, 7, 9, 15, 25)) +
  theme_blog() +
  labs(
    x = "K (matrix dimension)",
    y = "Divergent transitions (%)",
    title = "Divergence Rate by Matrix Dimension",
    subtitle = "One Bernoulli(0.5) zero pattern per K and 8,000 post-warmup transitions per fit"
  ) +
  theme(legend.position = "bottom")

Figure 1: Divergence rate by matrix dimension and LKJ concentration. Thin gold lines are the original sequential transforms and the red line is cvine_lookahead; missing points indicate initialization failures. The tri model has zero divergences in every panel. For qr2 at K = 15 and eta = 2, one chain has 100% divergences. The donut model remains below 1.5% throughout.

ggplot(sweep |> filter(!failed),
       aes(x = K, y = ess_per_1k_grad, group = model, color = family)) +
  geom_line(aes(linewidth = family, alpha = family, linetype = family)) +
  geom_point(size = 1.6) +
  facet_wrap(~eta_lab, nrow = 1, labeller = label_parsed) +
  scale_color_manual(values = fam_colors, name = NULL) +
  scale_linewidth_manual(guide = "none", values = c(
    "Forced zeros" = 0.6, "C-vine (look-ahead)" = 0.9, "Subspace (QR basis)" = 1.1,
    "Subspace (donut)" = 1.1, "Subspace (triangular basis)" = 1.3)) +
  scale_alpha_manual(guide = "none", values = c(
    "Forced zeros" = 0.65, "C-vine (look-ahead)" = 1, "Subspace (QR basis)" = 1,
    "Subspace (donut)" = 1, "Subspace (triangular basis)" = 1)) +
  scale_linetype_manual(guide = "none", values = c(
    "Forced zeros" = "solid", "C-vine (look-ahead)" = "solid", "Subspace (QR basis)" = "solid",
    "Subspace (donut)" = "22", "Subspace (triangular basis)" = "solid")) +
  scale_x_continuous(breaks = c(5, 7, 9, 15, 25)) +
  scale_y_log10() +
  theme_blog() +
  labs(
    x = "K (matrix dimension)",
    y = "min ESS per 1,000 gradients (log)",
    title = "Sampling Efficiency by Matrix Dimension",
    subtitle = "Higher is better; missing points indicate initialization failures"
  ) +
  theme(legend.position = "bottom")

Figure 2: Sampling efficiency (min bulk ESS over free correlations per 1,000 gradient evaluations, log scale) by dimension and concentration.

Outright failures: initialization could not find a feasible point after 100 attempts per chain.

model	K = 25	K = 15	K = 7	K = 9
carpenter	4 of 4 etas	—	—	—
cvine	4 of 4 etas	4 of 4 etas	—	—
cvine_lookahead	4 of 4 etas	—	—	—
ldl2	4 of 4 etas	—	—	—
ldl3	4 of 4 etas	—	—	—
pinkney	4 of 4 etas	—	—	—
schur	4 of 4 etas	4 of 4 etas	4 of 4 etas	4 of 4 etas

At \(\eta=1\), the entry-wise models have 15–18% divergences for \(K=5\) and 85–94% for \(K=7\) and \(K=9\). The uniform LKJ prior places more mass near the part of the parameter space where a solved Cholesky entry can exceed the remaining row norm. Four curves coincide in the \(\eta=1\) panel because carpenter, pinkney, ldl2, and ldl3 have the same density in their raw coordinates until floating-point differences separate their trajectories.

Initialization also becomes difficult for the sequential models. For these particular patterns, schur fails to initialize from \(K=7\) onward and cvine fails from \(K=15\) onward. The cvine_lookahead model initializes and samples through \(K=15\), but still fails to find a feasible default initialization at \(K=25\). The other entry-wise models also fail at \(K=25\). This depends on the zero pattern—schur initialized for every pattern in Section 5—but the problem becomes more likely as the number of constraints increases. Every initial value in the subspace models satisfies the zero constraints by construction.

The tri model has zero divergences in all 20 configurations. The donut model remains below 1.5% divergences throughout, with small counts for \(\eta\le2\) and none for \(\eta\ge3\). Through \(K=15\), the minimum bulk ESS for tri remains near the number of posterior draws. At \(K=25\), the minimum is taken over 157 free correlations and ranges from about 5,100 to 6,600 for tri. The donut model ranges from about 4,700 to 6,600. Both complete in seconds.

sub_tab <- sweep |>
  filter(model %in% c("tri", "donut")) |>
  mutate(cell = sprintf("%s (%s)%s",
                        format(round(min_ess), big.mark = ","),
                        round(ess_per_1k_grad),
                        ifelse(divergences > 0,
                               paste0(" [", divergences, " div]"), "")))

sub_wide <- sub_tab |>
  select(model, K, N_zero, eta, cell) |>
  pivot_wider(names_from = eta, values_from = cell,
              names_prefix = "eta = ") |>
  left_join(
    sub_tab |>
      group_by(model, K) |>
      summarize(`wall clock (s)` = sprintf("%.1f-%.1f",
                                           min(seconds), max(seconds)),
                .groups = "drop"),
    by = c("model", "K")
  ) |>
  mutate(`zero / free cells` = paste0(N_zero, " / ", choose(K, 2) - N_zero),
         .after = K) |>
  select(-N_zero) |>
  arrange(K, desc(model))

knitr::kable(
  sub_wide, align = c("l", "r", "c", rep("c", 4), "c"),
  caption = "tri and the donut (m_radius = 6) across the full sweep: min bulk ESS of 8,000 draws (ESS per 1,000 gradients in parentheses; divergent transitions in brackets where nonzero -- tri has none anywhere). Wall clock is the range across eta for 4 parallel chains."
)

tri and the donut (m_radius = 6) across the full sweep: min bulk ESS of 8,000 draws (ESS per 1,000 gradients in parentheses; divergent transitions in brackets where nonzero – tri has none anywhere). Wall clock is the range across eta for 4 parallel chains.

model	K	zero / free cells	eta = 1	eta = 2	eta = 3	eta = 4	wall clock (s)
tri	5	6 / 4	6,839 (127)	8,619 (183)	8,129 (192)	8,530 (204)	0.1-0.3
donut	5	6 / 4	2,364 (20) [23 div]	3,865 (66)	6,715 (120)	8,832 (172)	0.3-0.5
tri	7	7 / 14	6,922 (87)	7,572 (134)	10,480 (188)	10,832 (196)	0.4-0.5
donut	7	7 / 14	2,031 (19) [42 div]	5,796 (70)	7,128 (117)	8,430 (151)	0.5-0.8
tri	9	13 / 23	6,419 (87)	6,707 (118)	7,424 (132)	8,221 (147)	0.6-0.7
donut	9	13 / 23	5,315 (36) [54 div]	5,969 (74) [1 div]	6,977 (123)	6,834 (120)	0.7-1.5
tri	15	64 / 41	7,298 (67)	7,836 (134)	7,761 (127)	8,340 (149)	1.3-2.0
donut	15	64 / 41	5,046 (41) [99 div]	7,532 (79)	6,515 (82)	7,523 (122)	1.8-3.1
tri	25	143 / 157	5,081 (42)	5,837 (51)	5,742 (50)	6,553 (55)	7.0-8.9
donut	25	143 / 157	4,667 (29) [98 div]	5,547 (46) [2 div]	5,524 (46)	5,907 (50)	9.3-14.4

The QR basis remains feasible when the entry-wise models fail, but it becomes expensive as \(K\) grows. At \(K=25\), \(\eta=1\), qr2 takes 852 seconds compared with 8.9 seconds for tri. At \(K=15\), \(\eta=2\), one qr2 chain has 100% divergences and the minimum ESS is 7. The ldl3 model has a similar isolated failure at \(K=9\), \(\eta=2\), with a minimum ESS of 19. At \(K=25\), the minimum ESS for qr2 is about 3,500–4,100, compared with 5,100–6,600 for tri, and qr2 requires many more gradient evaluations. tri has the highest ESS per gradient in 14 of the 20 configurations, cvine_lookahead in four, and qr2 in two.

9. When to Use Which

Scenario	Recommendation
Structural zeros only	Use tri
Structural zeros at very large \(K\)	Use tri, but profile the per-row \(d\times d\) Cholesky factorizations
Zeros plus other known correlation values	Use the entry-wise method (Pinkney 2024a, 2024b); tri currently handles zeros only
Zeros plus bounds on other correlations	Use the entry-wise method, which already supports bounded transforms (Pinkney 2024a, 2024b)

For structural zeros alone, I do not see a reason to prefer one of the entry-wise parameterizations. In these benchmarks tri is more reliable and at least as efficient. The entry-wise method remains the established option here for more general constraints. The underlying look-ahead C-vine also accepts nonzero equalities and bounded free correlations, but this zeros-only adapter does not benchmark those cases.

10. Discussion

The entry-wise construction solves each zero constraint but leaves the boundary of the feasible row norm in the parameter space. This produces divergences for the smaller examples and initialization failures for the larger ones. cvine_lookahead narrows sequential partial-correlation intervals using exact \(3\times3\) feasibility projections and performs much better, but dense patterns can still leave a nonlocal constrained-completion boundary. The subspace construction removes that boundary because every parameter value maps to a row satisfying the constraints.

The choice of basis is important. A complement obtained directly from QR is valid, but its coordinates change with the earlier rows and sampling becomes less efficient. The triangular basis fixes the free coordinates to matrix entries and only solves for the constrained coordinates. In these benchmarks, that gives zero divergences without sacrificing effective sample size or run time.

There are two limitations I have not addressed here. The general look-ahead C-vine also supports nonzero equalities and bounded free correlations, but this bake-off’s interface specializes it to exact zeros and does not compare those targets. Also, the small Cholesky factorization performed for each constrained row should be profiled for matrices with dimensions in the hundreds.

Acknowledgments

The exact-constraint subspace idea, zeros as orthogonality constraints with the sample lifted from a lower-dimensional hemisphere through the orthogonal complement, including the \(\det R\) Jacobian factor, is due to Seth Axen (Axen 2023), whose radial prior (Axen 2022) the donut variant relies on. The tri model is that idea with a canonical triangular choice of complement basis and the determinant factored analytically. Thanks to Bob Carpenter for the compact zeros-only formulation benchmarked here as carpenter, and for the Discourse conversation that prompted this bake-off. The C-vine construction rests on the vine and partial-correlation results of Joe (2006) and Lewandowski et al. (2009).

Appendix: Stan Code for All Thirteen Parameterizations

tri

unitvec-stereo-corr-zeros.stan

functions {
  matrix cholesky_corr_constrain_stereo_jacobian(int K, vector raw,
                                              array[,] int zeros) {
    // Exact-constraint subspace geometry but with a canonical
    // triangular basis instead of Householder QR: basis vector k is a unit
    // at free column k with the row's forced columns filled in by
    // back-substitution through the zero constraints, then orthonormalized
    // by the Cholesky factor of its small Gram matrix. The basis is
    // anchored to matrix entries and varies smoothly with earlier rows.
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      int n_z = 0;
      while (zero_idx + n_z <= size(zeros)
             && zeros[zero_idx + n_z, 1] == i) {
        n_z += 1;
      }
      int d = i - 1 - n_z;

      // subsphere coordinates via the modified stereographic chart:
      // z = sqrt(t + 2) / (t + 1) * y with t = y'y, diagonal 1 / (t + 1)
      vector[d] z;
      real r = 1;
      if (d > 0) {
        vector[d] y = raw[raw_idx:(raw_idx + d - 1)];
        real t = dot_self(y);
        z = (sqrt(t + 2) / (t + 1)) * y;
        r = 1 / (t + 1);
        jacobian += 0.5 * (d - 2) * log(t + 2) - (d + 1) * log1p(t);
        raw_idx += d;
      }

      if (n_z == 0) {
        if (d > 0) {
          L[i, 1:d] = z';
        }
        L[i, i] = r;
      } else {
        array[n_z] int zc;
        for (m in 1:n_z) {
          zc[m] = zeros[zero_idx + m - 1, 2];
        }

        // the joint-delta term -0.5 * log det(Gram of zero-partner rows)
        // factors: splitting the partner matrix A by forced/free coords,
        // A_Z is triangular with the partners' diagonals, and by Sylvester
        // det(A'A) = prod(L[j,j])^2 * det(T'T)
        for (m in 1:n_z) {
          jacobian += -log(L[zc[m], zc[m]]);
        }

        if (d > 0) {
          array[d] int free_cols;
          {
            int fc = 0;
            int zp = 1;
            for (c in 1:(i - 1)) {
              if (zp <= n_z && zc[zp] == c) {
                zp += 1;
              } else {
                fc += 1;
                free_cols[fc] = c;
              }
            }
          }
          // triangular basis: unit at each free column, forced columns
          // filled by back-substitution so every column is orthogonal to
          // the zero-partner rows
          matrix[i - 1, d] T = rep_matrix(0, i - 1, d);
          for (k in 1:d) {
            int col = free_cols[k];
            T[col, k] = 1;
            for (m in 1:n_z) {
              int j = zc[m];
              if (j > col) {
                T[j, k] = -dot_product(L[j, 1:(j - 1)], T[1:(j - 1), k])
                          / L[j, j];
              }
            }
          }
          // T's free-column rows are the identity, so T'T = I + W'W with
          // W the forced-column rows only
          matrix[n_z, d] W;
          for (m in 1:n_z) {
            W[m, ] = T[zc[m], ];
          }
          // orthonormalize: v = T * S'^{-1} * z has ||v|| = ||z||;
          // diagonal(S) doubles as the det(T'T) half of the delta term
          matrix[d, d] S = cholesky_decompose(add_diag(crossprod(W), 1.0));
          jacobian += -sum(log(diagonal(S)));
          vector[d] t = mdivide_right_tri_low(z', S)';
          L[i, 1:(i - 1)] = (T * t)';
        }
        L[i, i] = r;
        zero_idx += n_z;
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_stereo_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

donut

unitvec-donut-corr-zeros.stan

functions {
  matrix cholesky_corr_constrain_donut_jacobian(int K, vector raw,
                                                array[,] int zeros,
                                                real m_radius) {
    // Exact-constraint subspace geometry with the triangular basis, but the
    // subsphere is parameterized by normalizing a free vector ("donut"):
    // one extra parameter per row, u = y / ||y||, with -0.5 * y'y keeping
    // the radius Gaussian and independent of the direction. The diagonal
    // is |u_last| (the two mirror hemispheres fold onto the same L), and
    // the measure conversion dz = |u_last| dsigma contributes log|u_last|
    // on top of the LKJ weight supplied by the model block.
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      int n_z = 0;
      while (zero_idx + n_z <= size(zeros)
             && zeros[zero_idx + n_z, 1] == i) {
        n_z += 1;
      }
      int d = i - 1 - n_z;

      vector[d + 1] y = raw[raw_idx:(raw_idx + d)];
      real s = norm2(y);
      vector[d + 1] u = y / s;
      // radial prior centered at m_radius repels the norm from the
      // origin singularity of the normalization (Axen, "A better unit
      // vector", Stan Discourse post 30); m_radius = 0 recovers the
      // implicit chi-distributed radius. The direction marginal, and
      // hence the Omega posterior, is unchanged.
      jacobian += -0.5 * square(s - m_radius) + log(abs(u[d + 1]));
      L[i, i] = abs(u[d + 1]);
      raw_idx += d + 1;

      if (n_z == 0) {
        if (d > 0) {
          L[i, 1:d] = u[1:d]';
        }
      } else {
        array[n_z] int zc;
        for (m in 1:n_z) {
          zc[m] = zeros[zero_idx + m - 1, 2];
        }

        // -0.5 * log det(Gram of zero-partner rows) from the joint delta
        matrix[i - 1, n_z] A;
        for (m in 1:n_z) {
          A[, m] = L[zc[m], 1:(i - 1)]';
        }
        jacobian += -sum(log(diagonal(cholesky_decompose(crossprod(A)))));

        if (d > 0) {
          array[d] int free_cols;
          {
            int fc = 0;
            int zp = 1;
            for (c in 1:(i - 1)) {
              if (zp <= n_z && zc[zp] == c) {
                zp += 1;
              } else {
                fc += 1;
                free_cols[fc] = c;
              }
            }
          }
          // triangular basis: unit at each free column, forced columns
          // filled by back-substitution so every column is orthogonal to
          // the zero-partner rows
          matrix[i - 1, d] T = rep_matrix(0, i - 1, d);
          for (k in 1:d) {
            int col = free_cols[k];
            T[col, k] = 1;
            for (m in 1:n_z) {
              int j = zc[m];
              if (j > col) {
                T[j, k] = -dot_product(L[j, 1:(j - 1)], T[1:(j - 1), k])
                          / L[j, j];
              }
            }
          }
          // orthonormalize: v = T * S'^{-1} * u[1:d] has ||v|| = ||u[1:d]||
          matrix[d, d] S = cholesky_decompose(crossprod(T));
          vector[d] t = mdivide_right_tri_low(u[1:d]', S)';
          L[i, 1:(i - 1)] = (T * t)';
        }
        zero_idx += n_z;
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
  real<lower=0> m_radius;                     // radial prior center; 0 = chi radius
}
parameters {
  // one extra parameter per row: dim d + 1 per row sums to
  // choose(K, 2) - N_zero + (K - 1)
  vector[choose(K, 2) - N_zero + K - 1] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_donut_jacobian(K, raw, zeros,
                                                                m_radius);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

carpenter

carpenter-corr-zeros.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }
  
  matrix cholesky_corr_zeros_jacobian(int D,
                                      vector raw,
                                      array[,] int zeros) {
    int raw_idx = 1;
    int zero_idx = 1;
    matrix[D, D] L = rep_matrix(0, D, D);
    for (i in 1:D) {
      real stick = 1;
      for (j in 1:(i - 1)) {
        if (is_zero(i, j, zeros, zero_idx)) {
          real b = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
          L[i, j] = -b / L[j, j];  // implies Omega[i, j] == 0
          zero_idx += 1;
        } else {
      real sqrt_stick = sqrt(stick);
          L[i, j] = lower_upper_bound_jacobian(raw[raw_idx], -sqrt_stick, sqrt_stick);  // inside stick
          raw_idx += 1;
          jacobian += 0;
        }
    stick -= L[i, j]^2;
      }
      L[i, i] = sqrt(stick);
    }
    return L;
  }

  real lkj_cholesky_corr_zeros_lpdf(matrix L,
                                    real nu,
                                    array[,] int zeros) {
    real lp = lkj_corr_cholesky_lpdf(L | nu);  // over-adjusts
    int N_zero = size(zeros);
    for (n in 1:N_zero) {                      
      int col_idx = zeros[n, 2];
      lp -= log(L[col_idx, col_idx]);          // correct over-adjustment
    }
    return lp;
  }
}
data {
  int<lower=2> D;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(D, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(D, 2) - N_zero] raw;           // raw parameters
}
transformed parameters {
  matrix[D, D] L_Omega = cholesky_corr_zeros_jacobian(D, raw, zeros);
}
model {
  L_Omega ~ lkj_cholesky_corr_zeros(eta, zeros);
}
generated quantities {
  matrix[D, D] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

pinkney

pinkney-corr-zeros2.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_outer_jacobian(int K, vector raw,
                                                array[,] int zeros) {
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      // remaining radius of row i: sqrt(1 - sum of squares so far)
      real l_ij_old = 1;
      for (j in 1:(i - 1)) {
        real l_ij_old_x_l_jj = l_ij_old * L[j, j];
        real b1 = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
        // how to derive the bounds
        // we know that the correlation value C is bound by
        // b1 - Ljj * Lij_old <= C <= b1 + Ljj * Lij_old
        // Now we want our bounds to be enforced too so
        // max(lb, b1 - Ljj * Lij_old) <= C <= min(ub, b1 + Ljj * Lij_old)
        // We have Lij_new = (C - b1) / Ljj
        real x;
        if (is_zero(i, j, zeros, zero_idx)) {
          x = -b1;  // Omega[i, j] = b1 + Lij * Ljj = 0
          zero_idx += 1;
        } else {
          real low = max({-l_ij_old_x_l_jj, -1 - b1});
          real up = min({l_ij_old_x_l_jj, 1 - b1});
          x = lower_upper_bound_jacobian(raw[raw_idx], low, up);
          raw_idx += 1;
        }
        L[i, j] = x / L[j, j];
        // free entries: change of variables x -> L[i, j]
        // zero entries: Carpenter's correction for the LKJ over-adjustment
        jacobian += -log(L[j, j]);
        l_ij_old *= sqrt(1 - (L[i, j] / l_ij_old)^2);
      }
      L[i, i] = l_ij_old;
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_outer_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

ldl

pinkney-corr-zeros-ldl.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_ldl_jacobian(int K, vector raw,
                                              array[,] int zeros) {
    // Omega = L * diag(D) * L' with unit lower-triangular L;
    // the Cholesky factor is diag_post_multiply(L, sqrt(D))
    matrix[K, K] L = diag_matrix(rep_vector(1, K));
    vector[K] D;
    vector[K] log_D;
    int raw_idx = 1;
    int zero_idx = 1;
    D[1] = 1;
    log_D[1] = 0;

    for (i in 2:K) {
      // log remaining radius^2 of row i: log(1 - sum_{k<j} L[i,k]^2 * D[k])
      real log_r = 0;
      for (j in 1:(i - 1)) {
        real b1 = dot_product(D[1:(j - 1)]' .* L[j, 1:(j - 1)],
                              L[i, 1:(j - 1)]);
        // Omega[i, j] = b1 + L[i, j] * D[j] and the correlation C satisfies
        // |C - b1| <= sqrt(D[j] * r), intersected with (-1, 1)
        real x;
        if (is_zero(i, j, zeros, zero_idx)) {
          x = -b1;  // Omega[i, j] = b1 + Lij * Dj = 0
          zero_idx += 1;
        } else {
          real b2 = exp(0.5 * (log_r + log_D[j]));
          real low = max({-b2, -1 - b1});
          real up = min({b2, 1 - b1});
          x = lower_upper_bound_jacobian(raw[raw_idx], low, up);
          raw_idx += 1;
        }
        L[i, j] = x / D[j];
        // free entries: change of variables x -> L_chol[i, j] = x / sqrt(D[j])
        // zero entries: Carpenter's correction -log(L_chol[j, j])
        jacobian += -0.5 * log_D[j];
        log_r = log_diff_exp(log_r, log_D[j] + 2 * log(abs(L[i, j])));
      }
      log_D[i] = log_r;
      D[i] = exp(log_r);
    }
    return diag_post_multiply(L, sqrt(D));
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_ldl_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

ldl2

pinkney-corr-zeros-ldl2.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_ldl2_jacobian(int K, vector raw,
                                               array[,] int zeros) {
    // Omega = L * diag(D) * L' with unit lower-triangular L;
    // the Cholesky factor is diag_post_multiply(L, sqrt(D)).
    // Same as ldl version except the remaining radius is updated with
    // log1m of the ratio instead of log_diff_exp.
    matrix[K, K] L = diag_matrix(rep_vector(1, K));
    vector[K] D;
    vector[K] log_D;
    int raw_idx = 1;
    int zero_idx = 1;
    D[1] = 1;
    log_D[1] = 0;

    for (i in 2:K) {
      // log remaining radius^2 of row i: log(1 - sum_{k<j} L[i,k]^2 * D[k])
      real log_r = 0;
      for (j in 1:(i - 1)) {
        real b1 = dot_product(D[1:(j - 1)]' .* L[j, 1:(j - 1)],
                              L[i, 1:(j - 1)]);
        // Omega[i, j] = b1 + L[i, j] * D[j] and the correlation C satisfies
        // |C - b1| <= sqrt(D[j] * r), intersected with (-1, 1)
        real x;
        if (is_zero(i, j, zeros, zero_idx)) {
          x = -b1;  // Omega[i, j] = b1 + Lij * Dj = 0
          zero_idx += 1;
        } else {
          real b2 = exp(0.5 * (log_r + log_D[j]));
          real low = max({-b2, -1 - b1});
          real up = min({b2, 1 - b1});
          x = lower_upper_bound_jacobian(raw[raw_idx], low, up);
          raw_idx += 1;
        }
        L[i, j] = x / D[j];
        // free entries: change of variables x -> L_chol[i, j] = x / sqrt(D[j])
        // zero entries: Carpenter's correction -log(L_chol[j, j])
        jacobian += -0.5 * log_D[j];
        log_r += log1m(D[j] * L[i, j]^2 / exp(log_r));
      }
      log_D[i] = log_r;
      D[i] = exp(log_r);
    }
    return diag_post_multiply(L, sqrt(D));
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_ldl2_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

ldl3

pinkney-corr-zeros-ldl3.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_ldl3_jacobian(int K, vector raw,
                                               array[,] int zeros) {
    // Omega = L * diag(D) * L' with unit lower-triangular L;
    // the Cholesky factor is diag_post_multiply(L, sqrt(D)).
    // Same as the ldl versions except the remaining radius is tracked on
    // the natural scale: r -= D[j] * L[i, j]^2.
    matrix[K, K] L = diag_matrix(rep_vector(1, K));
    vector[K] D;
    int raw_idx = 1;
    int zero_idx = 1;
    D[1] = 1;

    for (i in 2:K) {
      // remaining radius^2 of row i: 1 - sum_{k<j} L[i,k]^2 * D[k]
      real r = 1;
      for (j in 1:(i - 1)) {
        real b1 = dot_product(D[1:(j - 1)]' .* L[j, 1:(j - 1)],
                              L[i, 1:(j - 1)]);
        // Omega[i, j] = b1 + L[i, j] * D[j] and the correlation C satisfies
        // |C - b1| <= sqrt(D[j] * r), intersected with (-1, 1)
        real x;
        if (is_zero(i, j, zeros, zero_idx)) {
          x = -b1;  // Omega[i, j] = b1 + Lij * Dj = 0
          zero_idx += 1;
        } else {
          real b2 = sqrt(r * D[j]);
          real low = max({-b2, -1 - b1});
          real up = min({b2, 1 - b1});
          x = lower_upper_bound_jacobian(raw[raw_idx], low, up);
          raw_idx += 1;
        }
        L[i, j] = x / D[j];
        // free entries: change of variables x -> L_chol[i, j] = x / sqrt(D[j])
        // zero entries: Carpenter's correction -log(L_chol[j, j])
        jacobian += -0.5 * log(D[j]);
        r -= D[j] * L[i, j]^2;
      }
      D[i] = r;
    }
    return diag_post_multiply(L, sqrt(D));
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_ldl3_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

cvine

cvine-corr-zeros.stan

functions {
  // traversal here is tree-by-tree (column major), so scan instead of the
  // merge-join used by the row-major models
  int is_zero_scan(int row_idx, int col_idx, array[,] int zeros) {
    for (n in 1:size(zeros)) {
      if (zeros[n, 1] == row_idx && zeros[n, 2] == col_idx) {
        return 1;
      }
    }
    return 0;
  }

  matrix corr_constrain_cvine_jacobian(int K, vector raw,
                                       array[,] int zeros, real eta) {
    matrix[K, K] C = identity_matrix(K);
    matrix[K, K] P = identity_matrix(K);
    // LKJ in C-vine form: partial corr in tree i ~ Beta(a - i/2 adjustments);
    // constant base a with the -0.5 * i * log1m(p^2) correction below
    real a = eta + 0.5 * (K - 1);
    int raw_idx = 1;

    // tree 1: marginal correlations with variable 1
    for (k in 2:K) {
      if (is_zero_scan(k, 1, zeros)) {
        P[1, k] = 0;  // prior terms at p = 0 are constants
      } else {
        P[1, k] = lower_upper_bound_jacobian(raw[raw_idx], -1, 1);
        raw_idx += 1;
        jacobian += beta_lpdf(0.5 * (P[1, k] + 1) | a, a)
                    - 0.5 * log1m(P[1, k]^2);
      }
      C[1, k] = P[1, k];
      C[k, 1] = P[1, k];
    }

    for (i in 2:(K - 1)) {
      for (j in (i + 1):K) {
        vector[i - 1] b1;
        vector[i - 1] b2;
        int m = i;
        for (k in 1:(i - 1)) {
          m -= 1;
          b1[k] = P[m, i] * P[m, j];
          b2[k] = sqrt((1 - P[m, i]^2) * (1 - P[m, j]^2));
        }
        // unwinding the partial correlation recursion is affine:
        // C[i, j] = A + B * P[i, j]
        real A = 0;
        real B = 1;
        for (k in 1:(i - 1)) {
          A = b1[k] + A * b2[k];
          B *= b2[k];
        }
        real p_ij;
        if (is_zero_scan(j, i, zeros)) {
          p_ij = -A / B;         // C[i, j] = 0
          jacobian += -log(B);   // delta in Omega coords -> partial coords
        } else {
          p_ij = lower_upper_bound_jacobian(raw[raw_idx], -1, 1);
          raw_idx += 1;
        }
        // LKJ prior on the tree-i partial (evaluated at p* for zeros)
        jacobian += beta_lpdf(0.5 * (p_ij + 1) | a, a)
                    - 0.5 * i * log1m(p_ij^2);
        P[i, j] = p_ij;
        C[i, j] = A + B * p_ij;
        C[j, i] = C[i, j];
      }
    }
    return C;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular index pairs
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  // prior is included via the beta terms in the constraint, so the
  // model block stays empty
  matrix[K, K] Omega = corr_constrain_cvine_jacobian(K, raw, zeros, eta);
}
model {
}

cvine_lookahead

cvine-lookahead-corr-zeros.stan

/*
  Look-ahead C-vine parameterization of a correlation matrix with exact structural
  zeros. This is the equality-specialized interface to the look-ahead C-vine
  transform in https://gist.github.com/spinkney/7d994fb03079c9edb8a0e83a40d2faad.

  Each structural zero is an exact ordinary-correlation equality and consumes
  no element of `raw`. All other ordinary correlations have bounds (-1, 1).
  The transform accumulates the reduced raw-to-free-correlation Jacobian, and
  the model adds the LKJ determinant kernel with respect to those free cells.
*/

functions {
  // sqrt((1-a^2)(1-b^2)), evaluated stably on [-1,1]^2.
  real unit_circle_root_product(real a, real b) {
    real q = (1 - a) * (1 + a) * (1 - b) * (1 + b);
    return sqrt(fmax(0.0, q));
  }

  /*
    Exact projection of the 3 x 3 correlation elliptope. Let x be the
    correlation to project, with y in [y_lo, y_hi] and r in [r_lo, r_hi].
    The result is [minimum feasible x, maximum feasible x]. Degenerate point
    intervals are retained for exact zero constraints.
  */
  vector corr3_project_first(real y_lo_in, real y_hi_in,
                             real r_lo_in, real r_hi_in) {
    real y_lo = fmax(-1.0, fmin(1.0, y_lo_in));
    real y_hi = fmax(-1.0, fmin(1.0, y_hi_in));
    real r_lo = fmax(-1.0, fmin(1.0, r_lo_in));
    real r_hi = fmax(-1.0, fmin(1.0, r_hi_in));
    real x_lo;
    real x_hi;
    vector[2] out;

    if (y_lo > y_hi || r_lo > r_hi) {
      reject("corr3_project_first received a reversed interval: y=[",
             y_lo, ",", y_hi, "], r=[", r_lo, ",", r_hi, "].");
    }

    // Maximum x.
    if (fmax(y_lo, r_lo) <= fmin(y_hi, r_hi)) {
      x_hi = 1.0;
    } else if (y_lo > r_hi) {
      x_hi = y_lo * r_hi + unit_circle_root_product(y_lo, r_hi);
    } else {
      x_hi = y_hi * r_lo + unit_circle_root_product(y_hi, r_lo);
    }

    // Minimum x.
    if (fmax(y_lo, -r_hi) <= fmin(y_hi, -r_lo)) {
      x_lo = -1.0;
    } else if (y_lo > -r_lo) {
      x_lo = y_lo * r_lo - unit_circle_root_product(y_lo, r_lo);
    } else {
      x_lo = y_hi * r_hi - unit_circle_root_product(y_hi, r_hi);
    }

    out[1] = fmax(-1.0, fmin(1.0, x_lo));
    out[2] = fmax(-1.0, fmin(1.0, x_hi));
    return out;
  }

  /*
    Map unconstrained free coordinates to a box-constrained correlation
    matrix while solving exact ordinary-correlation equalities. The complete
    reduced raw -> free ordinary-correlation Jacobian is accumulated here.
  */
  tuple(matrix, real) cvine_bound_equality_jacobian(
      vector raw,
      matrix corr_lower,
      matrix corr_upper,
      real interval_tolerance,
      real validation_tolerance) {
    int K = rows(corr_lower);
    int M_free = num_elements(raw);
    int pos_free = 1;
    array[K, K] int pair_is_free = rep_array(0, K, K);
    matrix[K, K] lo_work = corr_lower;
    matrix[K, K] hi_work = corr_upper;
    matrix[K, K] P_work = identity_matrix(K);
    matrix[K, K] C_work = identity_matrix(K);
    real logJ = 0.0;
    real logdet = 0.0;

    if (cols(corr_lower) != K
        || rows(corr_upper) != K
        || cols(corr_upper) != K) {
      reject("cvine_bound_equality_transform requires equally sized square ",
             "bound matrices.");
    }

    // Classify coordinates and prepare the initial conditional bounds.
    for (i in 1:K) {
      lo_work[i, i] = 1.0;
      hi_work[i, i] = 1.0;

      if (i > 1) {
        for (j in 1:(i - 1)) {
          real width = corr_upper[i, j] - corr_lower[i, j];
          real lo_ij = corr_lower[i, j];
          real hi_ij = corr_upper[i, j];

          if (width > 0.0) {
            pair_is_free[i, j] = 1;
            pair_is_free[j, i] = 1;
            lo_ij = fmax(-1.0, corr_lower[i, j]);
            hi_ij = fmin(1.0, corr_upper[i, j]);
          } else if (width == 0.0) {
            // Preserve a point equality exactly.
            lo_ij = corr_lower[i, j];
            hi_ij = corr_lower[i, j];
          } else {
            reject("cvine_bound_equality_transform received a reversed ",
                   "interval at pair (", i, ",", j, "): [",
                   corr_lower[i, j], ",", corr_upper[i, j], "].");
          }

          lo_work[i, j] = lo_ij;
          lo_work[j, i] = lo_ij;
          hi_work[i, j] = hi_ij;
          hi_work[j, i] = hi_ij;
        }
      }
    }

    // Sequential C-vine construction.
    for (i in 1:(K - 1)) {
      for (j in (i + 1):K) {
        int is_free = pair_is_free[i, j];
        real lo_direct = lo_work[i, j];
        real hi_direct = hi_work[i, j];
        real lo = fmax(-1.0, lo_direct);
        real hi = fmin(1.0, hi_direct);
        real x;
        real log_scale = 0.0;

        // d C[i,j] / d P[i,j], holding earlier C-vine rows fixed.
        if (i > 1) {
          for (m in 1:(i - 1)) {
            log_scale += 0.5 * (
                log1m(square(P_work[m, i]))
                + log1m(square(P_work[m, j])));
          }
        }

        if (is_free == 1) {
          real width;

          // Intersect the direct interval with every 3 x 3 projection.
          if (j < K) {
            for (k in (j + 1):K) {
              vector[2] projected = corr3_project_first(
                  lo_work[i, k], hi_work[i, k],
                  lo_work[j, k], hi_work[j, k]);
              lo = fmax(lo, projected[1]);
              hi = fmin(hi, projected[2]);
            }
          }

          width = hi - lo;
          if (!(width > interval_tolerance)) {
            reject("Empty free C-vine interval at partial pair (", i, ",", j,
                   "). Interval after look-ahead=[", lo, ",", hi, "].");
          }
          if (pos_free > M_free) {
            reject("cvine_bound_equality_transform received too few raw ",
                   "coordinates; missing the coordinate for pair (", i, ",",
                   j, ").");
          }

          {
            real q = inv_logit(raw[pos_free]);

            if (!(q > 0.0 && q < 1.0)) {
              reject("inv_logit(raw[", pos_free,
                     "]) numerically reached an endpoint; raw=",
                     raw[pos_free], ".");
            }

            x = lo + width * q;
            logJ += log(width)
                    + log_inv_logit(raw[pos_free])
                    + log1m_inv_logit(raw[pos_free])
                    + log_scale;
          }

          pos_free += 1;
        } else {
          real fixed_x = 0.5 * (lo_direct + hi_direct);

          if (abs(hi_direct - lo_direct) > validation_tolerance) {
            reject("Internal equality interval acquired nonzero width at pair (",
                   i, ",", j, "): [", lo_direct, ",", hi_direct, "].");
          }
          if (!(fixed_x > -1.0 && fixed_x < 1.0)) {
            reject("Exact equality implies a non-interior partial correlation ",
                   "at pair (", i, ",", j, "): ", fixed_x, ".");
          }

          // Validate rather than alter the exact point during look-ahead.
          if (j < K) {
            for (k in (j + 1):K) {
              vector[2] projected = corr3_project_first(
                  lo_work[i, k], hi_work[i, k],
                  lo_work[j, k], hi_work[j, k]);

              if (fixed_x < projected[1] - validation_tolerance
                  || fixed_x > projected[2] + validation_tolerance) {
                reject("Exact equality partial at pair (", i, ",", j,
                       ") is incompatible with future pair ", k,
                       ". Fixed partial=", fixed_x,
                       ", projected interval=[", projected[1], ",",
                       projected[2], "].");
              }
            }
          }

          x = fixed_x;
        }

        P_work[i, j] = x;
        P_work[j, i] = x;
        logdet += log1m(square(x));

        // Update future intervals in the current C-vine row.
        if (j < K) {
          for (k in (j + 1):K) {
            vector[2] future_projected = corr3_project_first(
                x, x, lo_work[j, k], hi_work[j, k]);

            if (pair_is_free[i, k] == 1) {
              real next_lo = fmax(lo_work[i, k], future_projected[1]);
              real next_hi = fmin(hi_work[i, k], future_projected[2]);

              if (!(next_hi - next_lo > interval_tolerance)) {
                reject("Future free C-vine interval became empty after partial (",
                       i, ",", j, "). Future pair=", i, ",", k,
                       ", interval=[", next_lo, ",", next_hi, "].");
              }

              lo_work[i, k] = next_lo;
              lo_work[k, i] = next_lo;
              hi_work[i, k] = next_hi;
              hi_work[k, i] = next_hi;
            } else {
              real fixed_future = 0.5 * (
                  lo_work[i, k] + hi_work[i, k]);

              if (abs(hi_work[i, k] - lo_work[i, k])
                  > validation_tolerance) {
                reject("Internal future equality interval acquired width at ",
                       "pair (", i, ",", k, ").");
              }
              if (fixed_future < future_projected[1] - validation_tolerance
                  || fixed_future > future_projected[2]
                                      + validation_tolerance) {
                reject("Sampling partial (", i, ",", j,
                       ") made exact future equality pair (", i, ",", k,
                       ") infeasible. Fixed future partial=", fixed_future,
                       ", projected interval=[", future_projected[1], ",",
                       future_projected[2], "].");
              }

              lo_work[i, k] = fixed_future;
              lo_work[k, i] = fixed_future;
              hi_work[i, k] = fixed_future;
              hi_work[k, i] = fixed_future;
            }
          }
        }
      }

      // Propagate remaining intervals to the next C-vine level.
      if (i < K - 1) {
        for (j in (i + 1):(K - 1)) {
          for (k in (j + 1):K) {
            real pij = P_work[i, j];
            real pik = P_work[i, k];
            real denom = sqrt((1 - square(pij)) * (1 - square(pik)));

            if (!(denom > 0.0)) {
              reject("Zero conditional-correlation denominator at level ", i,
                     " for pair (", j, ",", k, ").");
            }

            if (pair_is_free[j, k] == 1) {
              real next_lo = fmax(
                  -1.0, (lo_work[j, k] - pij * pik) / denom);
              real next_hi = fmin(
                  1.0, (hi_work[j, k] - pij * pik) / denom);

              if (!(next_hi - next_lo > interval_tolerance)) {
                reject("Conditional free box became empty after C-vine row ",
                       i, " for pair (", j, ",", k, "). Interval=[",
                       next_lo, ",", next_hi, "].");
              }

              lo_work[j, k] = next_lo;
              lo_work[k, j] = next_lo;
              hi_work[j, k] = next_hi;
              hi_work[k, j] = next_hi;
            } else {
              real fixed_current = 0.5 * (
                  lo_work[j, k] + hi_work[j, k]);
              real fixed_next;

              if (abs(hi_work[j, k] - lo_work[j, k])
                  > validation_tolerance) {
                reject("Internal equality interval acquired width before ",
                       "conditioning at pair (", j, ",", k, ").");
              }

              fixed_next = (fixed_current - pij * pik) / denom;
              if (!(fixed_next > -1.0 && fixed_next < 1.0)) {
                reject("Exact ordinary equality implies a non-interior ",
                       "conditional correlation after level ", i,
                       " for pair (", j, ",", k, "): ", fixed_next, ".");
              }

              lo_work[j, k] = fixed_next;
              lo_work[k, j] = fixed_next;
              hi_work[j, k] = fixed_next;
              hi_work[k, j] = fixed_next;
            }
          }
        }
      }
    }

    if (pos_free != M_free + 1) {
      reject("cvine_bound_equality_transform received ", M_free,
             " raw coordinates but consumed ", pos_free - 1, ".");
    }

    // Convert C-vine partial correlations to ordinary correlations.
    for (i in 1:(K - 1)) {
      for (j in (i + 1):K) {
        real value = P_work[i, j];

        if (i > 1) {
          int m = i;
          for (step in 1:(i - 1)) {
            m -= 1;
            value = P_work[m, i] * P_work[m, j]
                    + value
                      * sqrt((1 - square(P_work[m, i]))
                             * (1 - square(P_work[m, j])));
          }
        }

        C_work[i, j] = value;
        C_work[j, i] = value;
      }
    }

    jacobian += logJ;
    return (C_work, logdet);
  }
}

data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // LKJ concentration
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular index pairs
}

transformed data {
  matrix[K, K] corr_lower = rep_matrix(-1.0, K, K);
  matrix[K, K] corr_upper = rep_matrix(1.0, K, K);

  for (i in 1:K) {
    corr_lower[i, i] = 1.0;
    corr_upper[i, i] = 1.0;
  }

  if (N_zero > 0) {
    for (n in 1:N_zero) {
      int i = zeros[n, 1];
      int j = zeros[n, 2];

      if (i < 2 || i > K || j < 1 || j >= i) {
        reject("zeros[", n, "] must be a lower-triangular index pair; got (",
               i, ",", j, ").");
      }
      if (n > 1) {
        for (m in 1:(n - 1)) {
          if (zeros[m, 1] == i && zeros[m, 2] == j) {
            reject("Duplicate structural-zero pair (", i, ",", j, ").");
          }
        }
      }

      corr_lower[i, j] = 0.0;
      corr_lower[j, i] = 0.0;
      corr_upper[i, j] = 0.0;
      corr_upper[j, i] = 0.0;
    }
  }
}

parameters {
  // Exact zeros are omitted from the parameter vector.
  vector[choose(K, 2) - N_zero] raw;
}

transformed parameters {
  corr_matrix[K] Omega;
  real log_det_Omega;

  {
    tuple(matrix[K, K], real) transform = cvine_bound_equality_jacobian(
        raw, corr_lower, corr_upper, 1e-12, 1e-10);
    Omega = transform.1;
    log_det_Omega = transform.2;
  }
}

model {
  // LKJ density with respect to the free ordinary-correlation coordinates.
  target += (eta - 1.0) * log_det_Omega;
}

rowscale

pinkney-corr-zeros-rowscale.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_rowscale_jacobian(int K, vector raw,
                                                   array[,] int zeros) {
    // Samples L[i, j] directly (bounds pre-divided by L[j, j]) and keeps the
    // remaining radius in the untouched tail of the row: each entry of
    // L[i, j:i] holds the current radius until column j is sampled, then the
    // tail is scaled down. L[i, i] ends up as the diagonal automatically.
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      L[i, 1:i] = rep_row_vector(1, i);
      for (j in 1:(i - 1)) {
        real l_ij_old = L[i, j];  // remaining radius of row i
        real b1 = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
        if (is_zero(i, j, zeros, zero_idx)) {
          L[i, j] = -b1 / L[j, j];  // Omega[i, j] = b1 + Lij * Ljj = 0
          jacobian += -log(L[j, j]);  // Carpenter's correction
          zero_idx += 1;
        } else {
          // |L[i, j]| <= radius and Omega[i, j] = b1 + Lij * Ljj in (-1, 1);
          // sampling L directly, the 1 / L[j, j] factor is inside the
          // bound width, so no separate jacobian term
          real low = max({-l_ij_old, (-1 - b1) / L[j, j]});
          real up = min({l_ij_old, (1 - b1) / L[j, j]});
          L[i, j] = lower_upper_bound_jacobian(raw[raw_idx], low, up);
          raw_idx += 1;
        }
        L[i, (j + 1):i] *= sqrt(1 - (L[i, j] / l_ij_old)^2);
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_rowscale_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

schur

pinkney-corr-zeros-schur.stan

functions {
  int is_zero(int row_idx, int col_idx, array[,] int zeros, int zero_idx) {
    return zero_idx <= size(zeros)
      && zeros[zero_idx, 1] == row_idx
      && zeros[zero_idx, 2] == col_idx;
  }

  matrix cholesky_corr_constrain_schur_jacobian(int K, vector raw,
                                                array[,] int zeros) {
    // Canonical partial correlation construction: each free entry is the
    // remaining radius times tanh(raw), so |L[i, j]| < radius automatically
    // and no correlation bounds are needed. The untouched tail of each row
    // holds the current radius (row-scaling trick).
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      L[i, 1:i] = rep_row_vector(1, i);
      for (j in 1:(i - 1)) {
        real r = L[i, j];  // remaining radius of row i
        if (is_zero(i, j, zeros, zero_idx)) {
          real b1 = dot_product(L[j, 1:(j - 1)], L[i, 1:(j - 1)]);
          L[i, j] = -b1 / L[j, j];  // Omega[i, j] = b1 + Lij * Ljj = 0
          jacobian += -log(L[j, j]);  // Carpenter's correction
          zero_idx += 1;
          L[i, (j + 1):i] *= sqrt(1 - (L[i, j] / r)^2);
        } else {
          real x = raw[raw_idx];
          real cosh_x = cosh(x);
          // tanh jacobian and x -> L[i, j] = r * tanh(x);
          // the radius update sqrt(1 - tanh^2) is 1 / cosh
          jacobian += -2 * log(cosh_x) + log(r);
          L[i, j] = r * tanh(x);
          raw_idx += 1;
          L[i, (j + 1):i] /= cosh_x;
        }
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_schur_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

qr

qr-corr-zeros.stan

functions {
  matrix cholesky_corr_constrain_unitvec_jacobian(int K, vector raw,
                                                  array[,] int zeros) {
    // Row i of a Cholesky correlation factor is a unit vector in R^i with
    // positive last coordinate. Omega[i, j] = 0 means row i is orthogonal
    // to row j, so row i lives on the unit sphere in the orthogonal
    // complement of its zero-partner rows: the constraint is satisfied
    // exactly by construction and no proposal can ever be infeasible.
    //
    // Chart per row: unit vector u = (C * y, 1)' / s with s = sqrt(1 + y'y),
    // C an orthonormal basis of the feasible subspace. Jacobian onto the
    // subsphere is -(d + 2) * log(s); the joint delta over the row's zeros
    // contributes -0.5 * log det(Gram) = -sum(log|diag(R)|) from the QR.
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      int n_z = 0;
      while (zero_idx + n_z <= size(zeros)
             && zeros[zero_idx + n_z, 1] == i) {
        n_z += 1;
      }
      int d = i - 1 - n_z;

      if (n_z == 0) {
        vector[d] y = raw[raw_idx:(raw_idx + d - 1)];
        real s = sqrt(1 + dot_self(y));
        L[i, 1:(i - 1)] = y' / s;
        L[i, i] = 1 / s;
        jacobian += -(d + 2) * log(s);
        raw_idx += d;
      } else {
        // columns are the zero-partner rows restricted to 1:(i - 1)
        matrix[i - 1, n_z] A;
        for (n in 1:n_z) {
          A[, n] = L[zeros[zero_idx + n - 1, 2], 1:(i - 1)]';
        }
        // full qr: the complement basis lives in Q's last d columns and
        // diagonal(R) supplies the Gram term from one decomposition
        matrix[i - 1, i - 1] Q;
        matrix[i - 1, n_z] R;
        (Q, R) = qr(A);
        jacobian += -sum(log(abs(diagonal(R))));
        if (d > 0) {
          vector[d] y = raw[raw_idx:(raw_idx + d - 1)];
          real s = sqrt(1 + dot_self(y));
          L[i, 1:(i - 1)] = (Q[, (n_z + 1):(i - 1)] * y)' / s;
          L[i, i] = 1 / s;
          jacobian += -(d + 2) * log(s);
          raw_idx += d;
        } else {
          L[i, i] = 1;  // row fully constrained: e_i
        }
        zero_idx += n_z;
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_unitvec_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

qr2

qr2-corr-zeros.stan

functions {
  matrix cholesky_corr_constrain_unitvec2_jacobian(int K, vector raw,
                                                   array[,] int zeros) {
    // Same exact-constraint geometry as the unitvec version: row i lives on
    // the unit sphere in the orthogonal complement of its zero-partner rows,
    // so every raw value maps to a valid constrained matrix. The chart on
    // the subsphere is the canonical-partial-correlation (tanh) stick
    // construction instead of the gnomonic map, to lighten the tails.
    matrix[K, K] L = rep_matrix(0, K, K);
    int raw_idx = 1;
    int zero_idx = 1;
    L[1, 1] = 1;

    for (i in 2:K) {
      int n_z = 0;
      while (zero_idx + n_z <= size(zeros)
             && zeros[zero_idx + n_z, 1] == i) {
        n_z += 1;
      }
      int d = i - 1 - n_z;

      // subspace coordinates via tanh sticks: w_k = tanh(x_k) * radius,
      // radius update sqrt(1 - tanh^2) = 1 / cosh, so the whole chart
      // jacobian collapses to accumulated powers of log(cosh)
      vector[d] w;
      real r = 1;
      for (k in 1:d) {
        real x = raw[raw_idx];
        real cosh_x = cosh(x);
        w[k] = r * tanh(x);
        jacobian += -(d - k + 2) * log(cosh_x);
        r /= cosh_x;
        raw_idx += 1;
      }

      if (n_z == 0) {
        if (d > 0) {
          L[i, 1:d] = w';
        }
        L[i, i] = r;
      } else {
        // columns are the zero-partner rows restricted to 1:(i - 1);
        // -0.5 * log det(Gram) from the joint delta over this row's zeros
        matrix[i - 1, n_z] A;
        for (n in 1:n_z) {
          A[, n] = L[zeros[zero_idx + n - 1, 2], 1:(i - 1)]';
        }
        // full qr: the complement basis lives in Q's last d columns and
        // diagonal(R) supplies the Gram term from one decomposition
        matrix[i - 1, i - 1] Q;
        matrix[i - 1, n_z] R;
        (Q, R) = qr(A);
        jacobian += -sum(log(abs(diagonal(R))));
        if (d > 0) {
          L[i, 1:(i - 1)] = (Q[, (n_z + 1):(i - 1)] * w)';
        }
        L[i, i] = r;
        zero_idx += n_z;
      }
    }
    return L;
  }
}
data {
  int<lower=2> K;                             // dimension of correlation matrix
  real<lower=0> eta;                          // concentration in LKJ Cholesky
  int<lower=0, upper=choose(K, 2)> N_zero;    // # structural zero correlations
  array[N_zero, 2] int zeros;                 // lower triangular, row major order
}
parameters {
  vector[choose(K, 2) - N_zero] raw;
}
transformed parameters {
  matrix[K, K] L_Omega = cholesky_corr_constrain_unitvec2_jacobian(K, raw, zeros);
}
model {
  L_Omega ~ lkj_corr_cholesky(eta);
}
generated quantities {
  matrix[K, K] Omega = multiply_lower_tri_self_transpose(L_Omega);
}

References

Axen, Seth. 2022. A Better Unit Vector. Stan Discourse thread, post 30. https://discourse.mc-stan.org/t/a-better-unit-vector/26989/30.

Axen, Seth. 2023. Working with Positive Definite Matrices with Structural Zeros. TuringLang/Turing.jl GitHub discussion #2067. https://github.com/TuringLang/Turing.jl/discussions/2067.

Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, et al. 2017. “Stan: A Probabilistic Programming Language.” Journal of Statistical Software 76 (1): 1–32. https://doi.org/10.18637/jss.v076.i01.

Joe, Harry. 2006. “Generating Random Correlation Matrices Based on Partial Correlations.” Journal of Multivariate Analysis 97 (10): 2177–89. https://doi.org/10.1016/j.jmva.2005.05.010.

Lewandowski, Daniel, Dorota Kurowicka, and Harry Joe. 2009. “Generating Random Correlation Matrices Based on Vines and Extended Onion Method.” Journal of Multivariate Analysis 100 (9): 1989–2001. https://doi.org/10.1016/j.jmva.2009.04.008.

Pinkney, Sean. 2024a. A Short Note on a Flexible Cholesky Parameterization of Correlation Matrices. https://arxiv.org/abs/2405.07286.

Pinkney, Sean. 2024b. Structured Correlation Matrices. StanCon 2024. https://discourse.mc-stan.org/t/new-users-guide-intro-sem-and-user-defined-constrained-parameters-chapters/41430/4.

Pinkney, Sean. 2025. Updated Cholesky Corr Parameterization Testing. Stan Discourse thread. https://discourse.mc-stan.org/t/updated-cholesky-corr-parameterization-testing/38827.

Pinkney, Sean. 2026. A Better-Er Unit Vector. Stan Discourse thread. https://discourse.mc-stan.org/t/a-better-er-unit-vector/40830.

Footnotes

1. @ethothepie.bsky.social corrected the original here where I incorrectly implied the zeros in the correlation matrix are conditionally unrelated when this is marginal

2. The pass-zero schur fit failed during initialization for seed 5, while the one-pass transform initialized all four chains. Its divergence and efficiency comparison therefore uses the four patterns that succeeded under both settings.