Error with Bayesian Gaussian Process "PosDefException: matrix is not positive definite; Cholesky factorization failed."

I got the error of “PosDefException” when trying Bayesian Gaussian Process in Turing.jl. The Bayesian framework looks as follows

@model function GP_μ(μ, D; cov_fn::Function=exp_cov_fn)
    n_station = length(μ)

    # weakly informative prior over mean field for the μ_i
    μ_mean ~ Normal(100, 25)
	
    # weakly informative priors over the kernel variance parameter
    # kernel variance parameter is > 0
    μ_kernelvariance ~ InverseGamma(10, 50) # 10 ± 5 roughly

    # weakly informative priors over kernel length scale ℓ
    min_dist = minimum(D) + 1E-3
    max_dist = maximum(D)
    μ_kernellength ~ Uniform(min_dist, max_dist)
	
    ########################
    # GAUSSIAN PROCESS MODEL
    ########################

    # the kernels are computed exactly
    μ_kernel = cov_fn(D, μ_kernelvariance_log, μ_kernellength)
	
    # process layer implements the Gaussian Process
    δ = 1E-6 # add an offset for numerical stability
    μ ~ MvNormal(μ_mean * ones(n_station), Matrix(Hermitian(μ_kernel + δ * I)))

end

function calc_μ_true(lon, lat)
    dist = Euclidean()(CENTER, [lon, lat])
    return 100 + 25 * sin(2 * π * dist / 5 + 3π / 2) + rand(Normal(0, 10))/5
end

n_station = 30 # number of synthetic stations
n_obs = 30 # number of obs per station

# get the true values
lons = rand(Uniform(LON_RANGE[1], LON_RANGE[2]), n_station)
lats = rand(Uniform(LAT_RANGE[1], LAT_RANGE[2]), n_station)
μ_true = calc_μ_true.(lons, lats)
distance = calc_dist(hcat(lats, lons))

model_μ = GP_μ(μ_true, distance)
chains = begin
	mcmc_sampler = DynamicNUTS()
	n_per_chain = 1000
	nchains = 2
	sample(model_μ, mcmc_sampler, MCMCThreads(), n_per_chain, nchains; drop_warmup=true)
end

I tried adding a noise term, using Hermitian(), Symmetric(), but still got the error of either “matrix is not positive definite” or “matrix is not Hermitian” for the covariance I suppose. Any idea how I might solve this?

What do you get when you just compute the cholesky factorization using cholesky()?

I got results like this

LinearAlgebra.Cholesky{Float64, Matrix{Float64}}([7.712549336531022 1.847900191359362 … 1.4450590383187025 0.251966480687946; 14.252021394844194 7.4879023865966055 … 0.0032314205728996553 0.5100249881619099; … ; 11.145089127233064 2.694521435254619 … 4.61337750842342 -1.0599196726807111e-5; 1.9433039134578745 4.2846262339608705 … 0.45873498713978184 4.399661715685264], 'U', 0)

I got the PosDefException error from bayesian sampling, but computing cholesky() or isposdef() on the MvNormal covariance (μ_kernel) seemed fine.