# 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:** 1\
**Showing post:** 11

<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:** [August 19, 2023, 11:01am UTC](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020/11 "2023-08-19T11:01:12Z")

</div>

This implementation _should_ work, but it fails with the same `PosDefException`.

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

Random.seed!(121)

# generate data
sigma = [1,2,3]
Omega = [1 0.3 0.2;
         0.3 1 0.1;
         0.2 0.1 1]

Sigma = Diagonal(sigma) * Omega * Diagonal(sigma)
N = 100
J = 3

y = rand(MvNormal(Sigma), N)'

@model function correlation_chol(J, N, y)
    sigma ~ filldist(truncated(Cauchy(0., 5.); lower = 0), J)
    F ~ LKJCholesky(J, 1.0)
    Σ_L = Diagonal(sigma .+ 1e-3) * F.L
    Sigma = PDMat(Cholesky(Σ_L))
    for i in 1:N
        y[i,:] ~ MvNormal(Sigma)
    end
    return Sigma
end

# same ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
chain2 = sample(correlation_chol(J, N, y), NUTS(), 1000)

```

Here’s the full stracktrace.

> **Stacktrace**
>
> ```julia
> ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
> Stacktrace:
> [1] checkpositivedefinite
> @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-beta2+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:67 [inlined]
> [2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::NoPivot; check::Bool)
> @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-beta2+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:269
> [3] cholesky! (repeats 2 times)
> @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-beta2+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:267 [inlined]
> [4] cholesky(A::Hermitian{Float64, Matrix{Float64}}, ::NoPivot; check::Bool)
> @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-beta2+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:401
> [5] cholesky (repeats 2 times)
> @ Bijectors ~/.julia/juliaup/julia-1.10.0-beta2+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/cholesky.jl:401 [inlined]
> [6] cholesky_lower(X::Matrix{Float64})
> @ Bijectors ~/.julia/packages/Bijectors/fD1YU/src/utils.jl:31
> [7] transform(b::Bijectors.VecCholeskyBijector, X::Matrix{Float64})
> @ Bijectors ~/.julia/packages/Bijectors/fD1YU/src/bijectors/corr.jl:220
> [8] with_logabsdet_jacobian
> @ DynamicPPL ~/.julia/packages/Bijectors/fD1YU/src/bijectors/corr.jl:214 [inlined]
> [9] with_logabsdet_jacobian_and_reconstruct(f::Bijectors.VecCholeskyBijector, dist::LKJCholesky{…}, x::SubArray{…})
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/abstract_varinfo.jl:565
> [10] _inner_transform!(vi::DynamicPPL.TypedVarInfo{…}, vn::AbstractPPL.VarName{…}, dist::LKJCholesky{…}, f::Bijectors.VecCholeskyBijector)
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:885
> [11] macro expansion
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:785 [inlined]
> [12] _link!(metadata::@NamedTuple{…}, vi::DynamicPPL.TypedVarInfo{…}, vns::@NamedTuple{…}, ::Val{…})
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:771
> [13] _link!(metadata::@NamedTuple{…}, vi::DynamicPPL.TypedVarInfo{…}, vns::@NamedTuple{…}, ::Val{…})
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:769 [inlined]
> [14] _link!(metadata::@NamedTuple{…}, vi::DynamicPPL.TypedVarInfo{…}, vns::@NamedTuple{…}, ::Val{…})
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:765 [inlined]
> [15] link!!
> @ Turing.Inference ~/.julia/packages/DynamicPPL/W8pRQ/src/varinfo.jl:726 [inlined]
> [16] link!!
> @ Turing.Inference ~/.julia/packages/DynamicPPL/W8pRQ/src/abstract_varinfo.jl:383 [inlined]
> [17] initialstep(rng::TaskLocalRNG, model::DynamicPPL.Model{…}, spl::DynamicPPL.Sampler{…}, vi::DynamicPPL.TypedVarInfo{…}; init_params::Nothing, nadapts::Int64, kwargs::@Kwargs{})
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/hmc.jl:141
> [18] step(rng::TaskLocalRNG, model::DynamicPPL.Model{…}, spl::DynamicPPL.Sampler{…}; resume_from::Nothing, init_params::Nothing, kwargs::@Kwargs{…})
> @ DynamicPPL ~/.julia/packages/DynamicPPL/W8pRQ/src/sampler.jl:111
> [19] step
> @ AbstractMCMC ~/.julia/packages/DynamicPPL/W8pRQ/src/sampler.jl:84 [inlined]
> [20] macro expansion
> @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:125 [inlined]
> [21] macro expansion
> @ AbstractMCMC ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
> [22] (::AbstractMCMC.var"#21#22"{…})()
> @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:12
> [23] with_logstate(f::Function, logstate::Any)
> @ Base.CoreLogging ./logging.jl:515
> [24] with_logger(f::Function, logger::LoggingExtras.TeeLogger{Tuple{…}})
> @ Base.CoreLogging ./logging.jl:627
> [25] with_progresslogger(f::Function, _module::Module, logger::Logging.ConsoleLogger)
> @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:36
> [26] macro expansion
> @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:11 [inlined]
> [27] mcmcsample(rng::TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::@Kwargs{…})
> @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:116
> [28] sample(rng::TaskLocalRNG, model::DynamicPPL.Model{…}, sampler::DynamicPPL.Sampler{…}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, nadapts::Int64, discard_adapt::Bool, discard_initial::Int64, kwargs::@Kwargs{})
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/hmc.jl:121
> [29] sample
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/hmc.jl:91 [inlined]
> [30] #sample#3
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/Inference.jl:196 [inlined]
> [31] sample
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/Inference.jl:189 [inlined]
> [32] #sample#2
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/Inference.jl:186 [inlined]
> [33] sample(model::DynamicPPL.Model{…}, alg::NUTS{…}, N::Int64)
> @ Turing.Inference ~/.julia/packages/Turing/aBDzv/src/inference/Inference.jl:180
> [34] top-level scope
> @ REPL[28]:1
> Some type information was truncated. Use `show(err)` to see complete types.
> 
> ```

Since `cholesky` is not called in the constraining transform or in the model, the error must be happening in the \_un\_constraining transform. But since `F` is a `Cholesky`, I don’t see why a `cholesky` transform is being called at all. @torfjelde is there a missing overload for `Cholesky` or `CholeskyVariate` somewhere in DynamicPPL?

---

_[View the full topic](https://discourse.julialang.org/t/singular-exception-with-lkjcholesky/85020)._