Hyper Priors in Turing

Hi,

I’m wondering if hyperpriors work on Turing. I’m playing around with Turing and trying to use it to estimate a Hidden Markov Model. In my HMM model, I’ve used a multivariate normal prior on parameters of my transition probability model with its mean and variance components respectively having a multivariate normal and inverse wishart prior (hyperpriors). With the hyper priors, I always get the error:
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

Its tough to explain the code and model and be concise so I’ve created a similar example on bayesian linear regression with similar prior and hyper priors to estimate the regression coefficients. This time however, I did not get the earlier error seen above but was instead returned a long error as follows:

ERROR: MethodError: no method matching PDMats.PDMat{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#24")){Turing.Core.VarReplay.VarInfo,Turing.Model{Tuple{:σ₂,:Sigma,:Mu,:beta_hat},Tuple{:y},getfield(Main, Symbol("###inner_function#31011#140")).......... and more...

The code for the bayesian regression is as follows:

####################################################################
# Data Generation
####################################################################
df = DataFrame(
    x1 = rand(Normal(), 100), 
    x2 = rand(Normal(), 100)
);
beta = [5, 4, 5];
df[:y] = beta[1] .+ beta[2]*df[:x1] + beta[3]*df[:x2] + rand(Normal(), 100);


####################################################################
# Bayesian Inference with Turing
####################################################################
@model BayesReg(y, x1, x2) = begin
    
    n_obs = length(y)

    σ₂ ~ InverseGamma(0.001, 0.001)

    Sigma   ~ InverseWishart(
        3 + 1, 
        Matrix{Float64}(I, 3, 3) 
    )
    Mu ~ MvNormal( 
        zeros(3),    
        Matrix{Float64}(I, 3, 3)   
    )

    beta_hat ~ MvNormal(Mu, Sigma)

    mu = beta_hat[1] .+ beta_hat[2]*x1 .+ beta_hat[3]*x2
    for i = 1:n_obs
        y[i] ~ Normal(mu[i], σ₂)
    end
end;


model = BayesReg(df[:y], df[:x1], df[:x2])
chain = sample(model, NUTS(1500, 200, 0.65));

Any ideas where I may have gone wrong or misunderstood something about the package…?

I’m guessing the first error happens because the covariance matrix of your hyperprior is not quite symmetric due to floating point rounding error:

using Distributions, Random, LinearAlgebra
seed!(1)
A = randn(3,3);
A = A'A
B = randn(3,3)
B = B'B
C = A' * B * A
# Result looks symmetric...
# 3×3 Array{Float64,2}:
#  17.0807   -5.12167  -11.5349
# -5.12167   4.1515     5.32059
# -11.5349    5.32059    9.48419
# ...but isn't actually:
issymetric(C)
# false
MvNormal(zeros(3), C)
# ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

The solution is to wrap C in Symmetric to tell MvNormal it can ignore the floating-point errors and assume the matrix is in fact symmetrical:

MvNormal(zeros(3), Symmetric(C)) # works

This is a confusion that I’ve run into before, as have other people. It would be great to get at least a more informative error message for it.

The second error is coming from the automatic differentiation done by the NUTS algorithm–if you try

chain = sample(model, MH(1500))

it runs. I can’t parse that error message either, but it looks like there’s some unhappy interaction going on between types from PDMats, ForwardDiff, and Turing.

3 Likes

Thanks @ElOceanografo, you are completely right about the issue with symmetry.

You are also completely right about the issue with the automatic differentiation (AD). I managed to resolve it by changing the AD to Turing.setadbackend(:reverse_diff).

Can’t thank you enough, your answer was prompt and came with examples. Was about to give up on Turing!

1 Like

No worries, glad you sorted it out!

Pinging @andreasnoack and @jrevels also…can you tell if the issue with ForwardDiff/PDMats here is something expected, or is it a bug?

I found the problem, but while fixing it, I may have stumbled upon a bigger problem in Julia itself; it silently exits. I will open an issue in Julia.

Edit: https://github.com/JuliaLang/julia/issues/31067

Hi all,

Sorry to chime in here, but I am having the same issue. However, I’ve found that Turing only throws ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed. when the dimensionality of the covariance matrix increases above 3. For exampe, I tried @krenova’s MWE (3 dimensional Sigma) and everything works. However, when I slightly modify the code to include a 4th variable I get the PosDefException. Below is the modified code, including the Symmetric(...) fix suggested by @ElOceanografo.

df = DataFrame(x1 = rand(Normal(), 100), x2 = rand(Normal(), 100), x3 = rand(Normal(), 100))
beta = [5, 4, 5, 6]
df[!, :y] = beta[1] .+ beta[2]*df[!, :x1] + beta[3]*df[!, :x2] + beta[4]*df[!, :x3] + rand(Normal(), 100)

@model BayesReg(y, x1, x2, x3) = begin
    
    n_obs = length(y)

    σ₂ ~ InverseGamma(0.001, 0.001)

    Sigma   ~ InverseWishart(
        4 + 1, 
        Matrix{Float64}(I, 4, 4) 
    )
    Mu ~ MvNormal( 
        zeros(4),    
        Matrix{Float64}(I, 4, 4)   
    )

    beta_hat ~ MvNormal(Mu, Symmetric(Sigma))

    mu = beta_hat[1] .+ beta_hat[2]*x1 .+ beta_hat[3]*x2 .+ beta_hat[4]*x3 
    for i = 1:n_obs
        y[i] ~ Normal(mu[i], σ₂)
    end
end

model = BayesReg(df[!, :y], df[!, :x1], df[!, :x2], df[!, :x3])
chain = sample(model, NUTS(0.65), 1000)

I’m using Turing v0.15.10. Does anybody know how to fix this? Any help would be greatly appreciated!

Ps.The issue mentioned by @ElOceanografo contains a suggestion that GitHub - timholy/PositiveFactorizations.jl: Positive-definite "approximations" to matrices can be used to help with the apparent numerical issues, but I’m not sure how to actually do that within Turing…