Hyper Priors in Turing

Specific Domains Statistics
1

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

2

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
3

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
4

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?

5

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: Julia 1.1 silently exits · Issue #31067 · JuliaLang/julia · GitHub

6

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…