Singular Exception with LKJCholesky

Specific Domains Probabilistic programming
11

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

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 
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?

1 Like