Remove manual Cholesky workaround when PyMC fixes LKJ bug

[thomaspinder]

Context

PyMC 5.28 + PyTensor 2.38 + NumPy 2.4 have a compatibility bug where pm.LKJCholeskyCov and pm.LKJCorr fail with:

ValueError: not enough values to unpack (expected 5, got 3)

in pytensor/tensor/einsum.py. The einsum_path returns 3-tuples but PyTensor expects 5-tuples.

Current workaround

src/litterman/spec.py uses a manual Cholesky parameterization instead:

sd = pm.HalfCauchy("sigma_sd", beta=2.5, shape=n_vars)
L = pt.zeros((n_vars, n_vars))
L = pt.set_subtensor(L[np.diag_indices(n_vars)], sd)
if n_tril > 0:
    tril_vals = pm.Normal("tril_offdiag", mu=0, sigma=0.5, shape=n_tril)
    # fill lower triangle manually

What needs to happen

• Monitor PyMC/PyTensor releases for a fix to the LKJ einsum bug
• Once fixed, replace manual Cholesky with pm.LKJCholeskyCov for cleaner parameterization
• Test that the LKJ version produces equivalent posteriors

Labels

tech-debt, upstream-dependency

[thomaspinder]

Status update (2026-07-29): upstream pymc-devs/pymc#7101 was closed on 2026-04-03, but the fix postdates the release we pin — verified empirically just now on the installed stack (pymc 5.28.0 / pytensor 2.38.0 / numpy 2.4.2): LKJCholeskyCov still fails with the einsum unpacking error (ValueError: not enough values to unpack (expected 5, got 3)). So this remains blocked on a PyMC version bump. Note for whoever picks it up after the bump: removing the manual Cholesky workaround is a breaking parameterisation change — posterior variable names (sigma_sd, tril_offdiag) are load-bearing in the convergence-report block mapping (#176) and asserted by the Student-t branch's tests (#168), so migrate only after the current PR fleet merges, in a dedicated PR with a posterior-compatibility note.

[thomaspinder]

The 2026-07-29 adversarial audit turned up a finding that changes the shape of this migration: #253 — Constant-volatility prior on Σ depends on endogenous column order (src/impulso/volatility.py:88, P2).

The manual Cholesky parameterisation implies Sigma[i,i] = sd[i]² (1 + Σ_{j<i} tril²), so the prior variance of variable i is inflated by exactly 1 + i·tril_sigma² purely by its position in the endogenous column ordering (Monte-Carlo confirms 1.00/1.25/1.50/1.75 at n=4, and 3.25x at n=10). That is a non-exchangeable prior on an order-free object, and LKJCholeskyCov is exchangeable in the variable ordering.

So the assumption embedded in this issue — that the workaround is prior-equivalent to LKJCholeskyCov and swapping back is a pure refactor — is disproved. The eventual migration must be treated as a prior CHANGE, not a refactor: posteriors will move (by an amount that grows with n and with the position of the variable), so it needs a changelog entry, a note in the docs, and regression baselines refreshed rather than asserted unchanged.

Two possible sequencings, both fine:

• Fix Constant-volatility prior on Σ depends on endogenous column order #253 first (reparameterise so implied marginal SDs are position-invariant), which makes this issue a genuine refactor afterwards.
• Or leave Constant-volatility prior on Σ depends on endogenous column order #253 open and fold it into this migration, accepting the prior change at the swap point.

---

Found by the 2026-07-29 adversarial audit (reports/audit-2026-07-29.md); every finding was independently verified with a reproduction before filing.