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

stan
correlation
parameterizations
Author

Sean Pinkney

Published

July 15, 2026

Modified

July 25, 2026

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 marginally1 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 twelve 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 solves for individual Cholesky entries so that the requested correlations are zero. 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

The first group uses 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 eight 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)
rowscale direct Cholesky entries, row-scaling stick
schur canonical partial correlations, \(L_{ij} = r \tanh(x)\)

These 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, each variant has divergences in about 1% of post-warmup transitions.

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 twelve 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 rowscale ldl3 schur qr qr2 tri donut
seed 1 34 53 53 26 41 44 40 40 0 0 0 0
seed 2 61 65 74 107 64 78 104 60 0 0 0 0
seed 3 77 44 67 57 69 44 68 59 0 0 0 0
seed 4 137 157 130 112 186 155 152 123 0 0 0 0
seed 5 110 94 114 54 141 114 102 29 0 0 0 0
total 419 413 438 356 501 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 rowscale ldl3 schur qr qr2 tri donut
seed 1 11639 10566 9979 10618 10219 11319 9840 11036 3995 3362 11238 10676
seed 2 10817 9866 10826 9772 10188 10153 10361 10040 4049 3728 11923 9045
seed 3 9344 11014 10244 10314 10946 9420 10625 10787 4334 4304 10655 9429
seed 4 8395 9159 8326 9257 8771 9463 8266 9458 3701 3815 9400 8161
seed 5 10381 10010 9268 8760 10301 9218 9896 4770 4267 4153 11845 8619
median 10381 10010 9979 9772 10219 9463 9896 10040 4049 3815 11238 9045

Efficiency (medians over the five patterns):

carpenter pinkney ldl ldl2 cvine rowscale ldl3 schur qr qr2 tri donut
min ESS / 1000 gradients 188.5 179.6 185.9 182.2 188.9 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 4.21 4.51 4.15 2.07 2.06 1.86 2.39

The tri model has zero divergences for all five patterns. It also has 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. 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 six entry-wise 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.

Code 
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)",
      .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,
  "Subspace (QR basis)" = blog_colors$purple,
  "Subspace (donut)" = blog_colors$blue,
  "Subspace (triangular basis)" = blog_colors$teal
)
Code 
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, "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, "Subspace (QR basis)" = 1,
    "Subspace (donut)" = 1, "Subspace (triangular basis)" = 1)) +
  scale_linetype_manual(guide = "none", values = c(
    "Forced zeros" = "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. Each thin gold line is one entry-wise parameterization; 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% at eta = 1 and has zero divergences otherwise.
Code 
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, "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, "Subspace (QR basis)" = 1,
    "Subspace (donut)" = 1, "Subspace (triangular basis)" = 1)) +
  scale_linetype_manual(guide = "none", values = c(
    "Forced zeros" = "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 — —
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 13–17% 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 entry-wise models. For these particular patterns, schur fails to initialize from \(K=7\) onward and cvine fails at \(K=15\). At \(K=25\), none of the six entry-wise models finds a feasible initial value in 100 attempts per chain. 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 has a small number at \(\eta=1\) and none for \(\eta\ge2\). 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.

Code 
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 17 of the 20 configurations; qr2 is higher in three of the \(K=5\) configurations.

7. 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 useful because it handles more general constraints: correlations fixed at nonzero values and bounds on otherwise free correlations.

8. 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. 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. Adding bounds on otherwise free correlations changes the row’s feasible region from a sphere to the intersection of a sphere and a box. 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 Twelve Parameterizations

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);
}
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-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-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);
}
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);
}
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);
}
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-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 {
}
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);
}
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-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-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↩︎