# Singular Exception with LKJCholesky

**URL:** <https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [July 30, 2022, 3:56am UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020 "2022-07-30T03:56:17Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [September 1, 2023, 7:00pm UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/21 "2023-09-01T19:00:52Z")

</div>

Alrighty! It took a bit longer than I originally thought because I realized we should clean up some stuff along the way, but here is the PR: [Remove usage of `DynamicPPL.tonamedtuple` by torfjelde · Pull Request #2071 · TuringLang/Turing.jl · GitHub](https://github.com/TuringLang/Turing.jl/pull/2071)

Once that’s through, `LKJCholeksy` (and others) should also work nicely with MCMChcains, etc.

---

<div class="post-metadata">

**Author:** ![KotaAihara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kotaaihara/32/43089_2.png) [@KotaAihara](https://discourse.julialang.org/u/KotaAihara)\
**Post date:** [October 12, 2023, 5:40am UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/22 "2023-10-12T05:40:19Z")

</div>

I made 2-dimentional State stace model . But the posterior of the coefficient matrix deviates from what I had inputted. Would you have any insight?

```julia
#seed
Random.seed!(1234)
T = 50

B = [0.9 0.1; 0.0 0.9]
Σ_X = [0.2 0.1; 0.1 0.2]
Σ_Y = [0.1 0.05; 0.05 0.1]
X = zeros(T, 2)
X[1, :] = [0, 0]
Y = zeros(T, 2)
for t in 2:T
    X[t, :] = B * X[t-1, :] + rand(MvNormal(Σ_X))
    Y[t, :] = X[t, :] + rand(MvNormal(Σ_Y))
end

@model function SSM_2D(y, ::Type{TV}=Matrix{Float64}) where {TV}
    T = size(y, 1)
    k = size(y, 2)
    #prior
    ###----
    F_x ~ LKJCholesky(k, 1.0)
    #σ_x : 2 ×1 vector
    σ_x ~ filldist(truncated(Cauchy(0., 5.); lower = 0), k)
    Σ_xL = Diagonal(σ_x .+ 1e-3) * F_x.L
    Σ_x = PDMat(Cholesky(Σ_xL + eps() * I))

    F_y ~ LKJCholesky(k, 1.0)
    σ_y ~ filldist(truncated(Cauchy(0., 5.); lower = 0), k)
    Σ_yL = Diagonal(σ_y .+ 1e-3) * F_y.L
    Σ_y = PDMat(Cholesky(Σ_yL + eps() * I))
    ###---- 
    B = TV(undef, k, k)
    for i in 1:k
        for j in 1:k
            B[i,j] ~ Normal(0, 1)
        end
    end
    x = TV(undef, T, k)
    for t in 2:T
        # Debugging statements to check the types.
        x[t, :] ~ MvNormal(B * x[t-1, :], Σ_x)
        y[t, :] ~ MvNormal(x[t, :], Σ_y)
    end
end

m = SSM_2D(Y)
chain = sample(m, NUTS(0.8), 1000);

```

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [October 12, 2023, 10:57am UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/23 "2023-10-12T10:57:25Z")

</div>

Can you say more about how it deviated from your expectations and what those expectations were?

---

<div class="post-metadata">

**Author:** ![KotaAihara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kotaaihara/32/43089_2.png) [@KotaAihara](https://discourse.julialang.org/u/KotaAihara)\
**Post date:** [October 12, 2023, 5:11pm UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/24 "2023-10-12T17:11:31Z")

</div>

Here is my MCMC implement result. means of parameters is not true parameter and  
std is not estimated

> parameters mean std mcse ess\_bulk ess\_tail rhat  
> F\_x.L[1,1] 1.0000 0.0000 NaN NaN NaN NaN ⋯  
> F\_x.L[2,1] 0.8824 0.0000 0.0000 NaN NaN NaN ⋯  
> F\_x.L[2,2] 0.4706 0.0000 0.0000 NaN NaN NaN ⋯  
> σ\_x[1] 0.1356 0.0000 0.0000 NaN NaN NaN ⋯  
> σ\_x[2] 3.7133 0.0000 0.0000 NaN NaN NaN ⋯  
> F\_y.L[1,1] 1.0000 0.0000 NaN NaN NaN NaN ⋯  
> F\_y.L[2,1] 0.9602 0.0000 0.0000 NaN NaN NaN ⋯  
> F\_y.L[2,2] 0.2793 0.0000 0.0000 NaN NaN NaN ⋯  
> σ\_y[1] 5.6576 0.0000 0.0000 NaN NaN NaN ⋯  
> σ\_y[2] 5.5809 0.0000 0.0000 NaN NaN NaN ⋯  
> B[1,1] 1.4806 0.0000 0.0000 NaN NaN NaN ⋯  
> B[1,2] -0.3448 0.0000 0.0000 NaN NaN NaN ⋯  
> B[2,1] -0.7867 0.0000 0.0000 NaN NaN NaN ⋯  
> B[2,2] 1.8545 0.0000 0.0000 NaN NaN NaN ⋯

---

<div class="post-metadata">

**Author:** ![joshualeond](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joshualeond/32/230_2.png) [@joshualeond](https://discourse.julialang.org/u/joshualeond)\
**Post date:** [April 25, 2025, 4:41pm UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/25 "2025-04-25T16:41:08Z")

</div>

I know it has been years but just revisiting this today and was able to use your suggestion here and everything works well with this toy example. Here it is in case anyone is interested:

```julia
using LinearAlgebra, PDMats, Random, Turing

Random.seed!(121)

# generate data  
σ = [1,2,3]  
Ω = [1 0.3 0.2;  
         0.3 1 0.1;  
         0.2 0.1 1]

Σ = Diagonal(σ) * Ω * Diagonal(σ)  
N = 100  
J = 3
y = rand(MvNormal(Σ), N)'

@model function correlation_chol(J, N, y)
	σ ~ filldist(truncated(Cauchy(0., 5.); lower = 0), J)
	F ~ LKJCholesky(J, 1.0)
    Σ_L = Diagonal(σ) * F.L
    Σ = PDMat(Cholesky(Σ_L + eps() * I))
    for i in 1:N
      y[i,:] ~ MvNormal(Σ)
    end
    # Compute Omega
    Ω = F.L * F.L'
    return (Σ=Σ, Ω=Ω)
end

model = correlation_chol(J, N, y)
chain = sample(model, NUTS(), 1000)

# check Omega results
gen_quants = generated_quantities(model, chain)
omega_matrices = [gq.Ω for gq in gen_quants]
mean(omega_matrices)

```

Resulting in:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/a/daf9a2c56d24946acb688738758e52c00d052352.png)

and the following for Omega:

![image](https://global.discourse-cdn.com/julialang/original/3X/5/6/56ee1187ae56258021892d4ef4aef92245f3d2c8.png)

What’s inferred isn’t _perfect_ but I think on-par with the Stan example I originally referenced. Of course increasing `y` length increases the accuracy in the results as well.

Thanks everyone for all the work you do for the open source community. Excited to start playing with Turing again.

[Previous page](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020.md?page=1)