PosDef issues with InverseWishart?

[bratslavia]

I had been wanting to use the LKJ distribution, and ran into #144. However, when I decided to fall back to the usual Inverse Wishart, I seemed to have similar issues.

When I try to draw any samples whatsoever from an Inverse Wishart prior, I seem to get a PosDef error in Bijectors.jl - see MWE below:

using Turing, Distributions

@model cov_matrix(data) = begin
    N, P = size(data)
    Σ ~ InverseWishart(P, diagm(ones(P)))
end

m1 = sample(cov_matrix(rand(900, 10)), HMC(0.01, 20), 200)

Is there something I'm doing wrong? When I draw random samples from an Inverse Wishart manually using rand() and check for PD-ness, the samples are all fine. Versions for Distributions, Bijectors, and Turing are 0.24.15, 0.8.14, and 0.15.12, respectively.

Edit: Running a few more variations, it seems it fails under ForwardDiff using the code above (with 10 dimensions). If I use ReverseDiff, I get numerical errors during sampling (with basically all proposals rejected). And if I use a larger number of dimensions (e.g., using say 50 dimensions instead of 10), sampling with ReverseDiff fails a different way.

[bratslavia]

With a bit more testing, I don't really see reliable (frequent) errors using `reversediff (although as the dimension gets larger, the likelihood that it crashes with a posdef error during sampling is not super-rare, either).

With forwarddiff, it runs fine if the dimension is 3 or less. But as soon as the dimension hits 4 or more, every attempt to run the model above gives nothing but posdef failures:

``julia ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed. Stacktrace: [1] checkpositivedefinite(info::Int64) @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\factorization.jl:18 [2] cholesky!(A::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:Σ,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Σ, Tuple{}}, Int64}, Vector{InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}}, Vector{AbstractPPL.VarName{:Σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#13#14", (:data,), (), (), Tuple{Matrix{Float64}}, Tuple{}}, DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 8}}, ::Val{false}; check::Bool) @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\cholesky.jl:282 [3] #cholesky#130 @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\cholesky.jl:378 [inlined] [4] cholesky (repeats 2 times) @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\cholesky.jl:378 [inlined] [5] pd_logpdf_with_trans(d::InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}, X::Matrix{ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#1"{DynamicPPL.TypedVarInfo{NamedTuple{(:Σ,), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:Σ, Tuple{}}, Int64}, Vector{InverseWishart{Float64, PDMats.PDMat{Float64, Matrix{Float64}}}}, Vector{AbstractPPL.VarName{:Σ, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#13#14", (:data,), (), (), Tuple{Matrix{Float64}}, Tuple{}}, DynamicPPL.Sampler{HMC{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.UnitEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 8}}, transform::Bool) @ Bijectors ~.julia\packages\Bijectors\nVGqm\src\Bijectors.jl:222 [6] logpdf_with_trans @ ~.julia\packages\Bijectors\nVGqm\src\Bijectors.jl:125 [inlined]

[bonStats]

I've just encountered this exact issue for dimension 4

[anfesanz]

The same error for 4 or more dimensions