Problems with Cholesky in MvNormal using a g prior

[gregorsteiner]

Hello,

I am trying to implement a simple Bayesian Model Averaging procedure. When I use a g prior on the coefficient vector in a linear model, I keep getting this error: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

When I wrap the covariance matrix with Symmetric() or Hermitian(), I get this error instead: PosDefException: matrix is not positive definite; Cholesky factorization failed. However, when manually checking the eigenvalues, they are all positive.

The solutions suggested in here or here did not help. A similar problem has been solved here, but I am not sure how to use that solution within a Turing model.

My code is below, any help would be highly appreciated. Thank you!

using Random, Turing, LinearAlgebra
using MCMCChains, Plots, StatsPlots
using SpecialFunctions, Optim

# simulate data
n = 50
p = 10
beta = zeros(p)
beta[1:5] .= 1

Random.seed!(42)
X = rand(MvNormal(zeros(p), 10 * I), n)'
y = 10 .+ X * beta + rand(Normal(0, 2), n)

# create model
@model function linmod(y, X)
    p = size(X, 2)
    n = size(X, 1)
    g = min(n, p^2)

    # priors 
    σ² ~ Exponential(10)
    α ~ Normal(0, 100)

    Σ = σ² * g * inv(X'X)
    β ~ MvNormal(zeros(p), Σ)

    # model
    y ~ MvNormal(α .+ X*β, σ² * I)

end

model = linmod(y, X)
optimize(model, MAP())

[torfjelde]

Sorry, I just saw this. Did you manage to address the issue?

[gregorsteiner]

No worries, I managed to fix the issue in this case by wrapping the covariance matrix in PDMat(Symmetric( )) , but I have encountered a few versions of this problem. I should have probably posted this in Distributions.jl or StaticArrays.jl, but it seems that they are already aware of this issue: https://github.com/JuliaStats/Distributions.jl/issues/1826 and https://github.com/JuliaArrays/StaticArrays.jl/issues/1218. Many thanks!