Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PDMat on vector from arraydist fails with reversediff #2101

Open
tiemvanderdeure opened this issue Oct 13, 2023 · 3 comments
Open

PDMat on vector from arraydist fails with reversediff #2101

tiemvanderdeure opened this issue Oct 13, 2023 · 3 comments

Comments

@tiemvanderdeure
Copy link

I bumped into some more shenanigans with reversediff and covariance matrices. Generating a vector of numbers using arraydist, multiplying that with a matrix from LKJCholesky, and passing it back to PDMat to generate a covariance matrix fails.

The following model works with forwarddiff but errors with reversediff:

@model function model_with_arraydist()
    stds ~ arraydist([truncated(Normal(0, x); lower = 0.0) for x in [1,2,3]])
    F ~ LKJCholesky(3, 3.0)
    Sigma = PDMat(Cholesky(LowerTriangular(stds .* F.L)))
end

sample(model_with_arraydist(), HMC(0.1, 10), 50)

The same code with filldist instead of arraydist does not error

@model function model_with_filldist()
    stds ~ filldist(truncated(Normal(0, 0.5); lower = 0.0), 3)
    F ~ LKJCholesky(3, 3.0)
    Sigma = PDMat(Cholesky(LowerTriangular(stds .* F.L)))
end

sample(my_model(), HMC(0.1, 10), 50)

Stacktrace:

ERROR: MethodError: Cannot `convert` an object of type Matrix{ReverseDiff.TrackedReal{Float64, Float64, Nothing}} to an object of type ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}

Closest candidates are:
  convert(::Type{T}, ::T) where T<:ReverseDiff.TrackedArray
   @ ReverseDiff C:\Users\tsh371\.julia\packages\ReverseDiff\UJhiD\src\tracked.jl:270
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra C:\Users\tsh371\.julia\juliaup\julia-1.9.3+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
  convert(::Type{T}, ::T) where T<:AbstractArray
   @ Base abstractarray.jl:16
  ...

  [1] PDMat(mat::Matrix{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}, chol::Cholesky{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}})
    @ PDMats C:\Users\tsh371\.julia\packages\PDMats\VftST\src\pdmat.jl:16
  [2] PDMat(fac::Cholesky{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}})
    @ PDMats C:\Users\tsh371\.julia\packages\PDMats\VftST\src\pdmat.jl:20
  [3] macro expansion
    @ C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\compiler.jl:555 [inlined]
  [4] model_with_arraydist(__model__::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, __varinfo__::DynamicPPL.ThreadSafeVarInfo{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, 
ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{Set{DynamicPPL.Selector}}}}}, ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}, Vector{Base.RefValue{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}}}}, __context__::DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG})
    @ Main
  [5] _evaluate!!
    @ C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\model.jl:963 [inlined]
  [6] evaluate_threadsafe!!
    @ C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\model.jl:952 [inlined]
  [7] evaluate!!
    @ C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\model.jl:887 [inlined]
  [8] logdensity(f::LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, θ::ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}})
    @ DynamicPPL C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\logdensityfunction.jl:94
  [9] Fix1
    @ .\operators.jl:1108 [inlined]
 [10] ReverseDiff.GradientTape(f::Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}}, input::Vector{Float64}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
    @ ReverseDiff C:\Users\tsh371\.julia\packages\ReverseDiff\UJhiD\src\api\tape.jl:199
 [11] gradient!(result::DiffResults.MutableDiffResult{1, Float64, Tuple{Vector{Float64}}}, f::Function, input::Vector{Float64}, cfg::ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}})
    @ ReverseDiff C:\Users\tsh371\.julia\packages\ReverseDiff\UJhiD\src\api\gradients.jl:41
 [12] gradient!
    @ C:\Users\tsh371\.julia\packages\ReverseDiff\UJhiD\src\api\gradients.jl:41 [inlined]
 [13] logdensity_and_gradient(∇ℓ::LogDensityProblemsADReverseDiffExt.ReverseDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, Nothing}, x::Vector{Float64})
    @ LogDensityProblemsADReverseDiffExt C:\Users\tsh371\.julia\packages\LogDensityProblemsAD\JoNjv\ext\LogDensityProblemsADReverseDiffExt.jl:72
 [14] ∂logπ∂θ
    @ C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\hmc.jl:160 [inlined]
 [15] ∂H∂θ
    @ C:\Users\tsh371\.julia\packages\AdvancedHMC\dgxuI\src\hamiltonian.jl:38 [inlined]
 [16] phasepoint(h::AdvancedHMC.Hamiltonian{AdvancedHMC.UnitEuclideanMetric{Float64, Tuple{Int64}}, AdvancedHMC.GaussianKinetic, Base.Fix1{typeof(LogDensityProblems.logdensity), LogDensityProblemsADReverseDiffExt.ReverseDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, Nothing}}, Turing.Inference.var"#∂logπ∂θ#36"{LogDensityProblemsADReverseDiffExt.ReverseDiffLogDensity{LogDensityFunction{DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, DynamicPPL.SamplingContext{DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext, TaskLocalRNG}}, Nothing}}}, θ::Vector{Float64}, r::Vector{Float64})
    @ AdvancedHMC C:\Users\tsh371\.julia\packages\AdvancedHMC\dgxuI\src\hamiltonian.jl:80
 [17] phasepoint
    @ C:\Users\tsh371\.julia\packages\AdvancedHMC\dgxuI\src\hamiltonian.jl:159 [inlined]
 [18] initialstep(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, vi::DynamicPPL.TypedVarInfo{NamedTuple{(:stds, :F), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:stds, Setfield.IdentityLens}, Int64}, Vector{Product{Continuous, Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}, Vector{Truncated{Normal{Float64}, Continuous, Float64, Float64, Nothing}}}}, Vector{AbstractPPL.VarName{:stds, Setfield.IdentityLens}}, 
Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:F, Setfield.IdentityLens}, Int64}, Vector{LKJCholesky{Float64}}, Vector{AbstractPPL.VarName{:F, Setfield.IdentityLens}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}; init_params::Nothing, nadapts::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\hmc.jl:164
 [19] step(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, 
(), AdvancedHMC.UnitEuclideanMetric}}; resume_from::Nothing, init_params::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DynamicPPL C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\sampler.jl:111
 [20] step
    @ C:\Users\tsh371\.julia\packages\DynamicPPL\m0PXI\src\sampler.jl:84 [inlined]
 [21] macro expansion
    @ C:\Users\tsh371\.julia\packages\AbstractMCMC\fWWW0\src\sample.jl:125 [inlined]
 [22] macro expansion
    @ C:\Users\tsh371\.julia\packages\ProgressLogging\6KXlp\src\ProgressLogging.jl:328 [inlined]
 [23] macro expansion
    @ C:\Users\tsh371\.julia\packages\AbstractMCMC\fWWW0\src\logging.jl:9 [inlined]
 [24] mcmcsample(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ AbstractMCMC C:\Users\tsh371\.julia\packages\AbstractMCMC\fWWW0\src\sample.jl:116
 [25] sample(rng::TaskLocalRNG, model::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Turing.Inference C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:208
 [26] sample
    @ C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:197 [inlined]
 [27] #sample#2
    @ C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:194 [inlined]
 [28] sample
    @ C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:187 [inlined]
 [29] #sample#1
    @ C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:184 [inlined]
 [30] sample(model::DynamicPPL.Model{typeof(model_with_arraydist), (), (), (), Tuple{}, Tuple{}, DynamicPPL.DefaultContext}, alg::HMC{Turing.Essential.ReverseDiffAD{false}, (), AdvancedHMC.UnitEuclideanMetric}, N::Int64)
    @ Turing.Inference C:\Users\tsh371\.julia\packages\Turing\SsIU6\src\mcmc\Inference.jl:178
 [31] top-level scope
@torfjelde
Copy link
Member

This is unfortunately an issue with ReverseDiff.jl that is non-trivial to resolve. It comes from the fact that ReverseDiff.jl does not work well with "structural" arrays, e.g. PDMat or Cholesky.

It might be possible to re-use some internal methods we are using to avoid these issues, e.g. https://github.com/TuringLang/Bijectors.jl/blob/04b79dd46eca8cea2f988348c47bd5e720a2b9a4/ext/BijectorsReverseDiffExt.jl#L222-L230

That is, the following works on my end:

julia> using Turing, LinearAlgebra, ReverseDiff

julia> Turing.setadbackend(:reversediff);

julia> Turing.setrdcache(true);

julia> @model function model_with_arraydist()
           stds ~ arraydist([truncated(Normal(0, x); lower = 0.0) for x in [1,2,3]])
           F ~ LKJCholesky(3, 3.0)
           Sigma = Bijectors.pd_from_lower(stds .* F.L)
       end
model_with_arraydist (generic function with 2 methods)

julia> mean(sample(model_with_arraydist(), NUTS(0.65), 1000))
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling 100%|████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:34
Mean
  parameters      mean 
      Symbol   Float64 

     stds[1]    0.9353
     stds[2]    1.2679
     stds[3]    3.4438
    F.L[1,1]    1.0000
    F.L[2,1]   -0.1247
    F.L[3,1]    0.1726
    F.L[2,2]    0.9243
    F.L[3,2]    0.0878
    F.L[3,3]    0.9025


julia> @model function model_with_filldist()
           stds ~ filldist(truncated(Normal(0, 0.5); lower = 0.0), 3)
           F ~ LKJCholesky(3, 3.0)
           Sigma = Bijectors.pd_from_lower(stds .* F.L)
       end
model_with_filldist (generic function with 2 methods)

julia> mean(sample(model_with_filldist(), NUTS(0.65), 1000))
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling 100%|████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:30
Mean
  parameters      mean 
      Symbol   Float64 

     stds[1]    0.3786
     stds[2]    0.2528
     stds[3]    0.3489
    F.L[1,1]    1.0000
    F.L[2,1]    0.1822
    F.L[3,1]   -0.0433
    F.L[2,2]    0.9302
    F.L[3,2]    0.1201
    F.L[3,3]    0.8273

Bijectors.pd_from_lower is not documented as it's not meant to be a public-facing method, but it effectively takes in a Matrix, assumes it's lower triangular, and the constructs a PDMat from this, "hiding" the PDMat(Cholesky(...)) from ReverseDiff.jl.

Also, just for the record, both of the examples you give fail on my end (which is what I expected) 😕

@tiemvanderdeure
Copy link
Author

tiemvanderdeure commented Oct 16, 2023

I'll try with Bijectors.pd_from_lower, thanks!

Also, just for the record, both of the examples you give fail on my end (which is what I expected) 😕

I made a typo in my MWE, but the following code definitely runs for me (I'm on Turing v0.92.2 and ReverseDiff v1.15.1)

using Turing, PDMats, LinearAlgebra, ReverseDiff

Turing.setadbackend(:reversediff)

@model function model_with_filldist()
    stds ~ filldist(truncated(Normal(0, 0.5); lower = 0.0), 3)
    F ~ LKJCholesky(3, 3.0)
    Sigma = PDMat(Cholesky(LowerTriangular(stds .* F.L)))
end

sample(model_with_filldist(), HMC(0.1, 10), 50)

@tiemvanderdeure
Copy link
Author

Using Bijectors.pd_from_lower I quickly ran into numerical stability problems, especially with bigger covariance matrices. E.g.

@model function model_with_filldist(i)
    stds ~ filldist(truncated(Normal(0, 0.5); lower = 0.0), i)
    F ~ LKJCholesky(i, 3.0)
    Sigma = Bijectors.pd_from_lower(stds .* F.L)
    Y ~ MvNormal(Sigma)
end

mean(sample(model_with_filldist(3), NUTS(0.65), 1000)) # Works
mean(sample(model_with_filldist(5), NUTS(0.65), 1000)) # Errors

Where the error i get is PosDefException: matrix is not Hermitian; Cholesky factorization failed.

My quick fix is to wrap the Sigma in a Hermitian, which runs (but probably isn't great for ReverseDiff either?)

@model function model_with_filldist2(i)
    stds ~ filldist(truncated(Normal(0, 0.5); lower = 0.0), i)
    F ~ LKJCholesky(i, 3.0)
    Sigma = Hermitian(Bijectors.pd_from_lower(stds .* F.L))
    Y ~ MvNormal(Sigma)
end

mean(sample(model_with_filldist2(3), NUTS(0.65), 1000)) # Works
mean(sample(model_with_filldist2(5), NUTS(0.65), 1000)) # Works

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants