Hi All,
I’m using julia to sample from a gaussian process with squared exponential kernel.
using Distributions, Distances
x_star = -5:0.1:5
K = exp(-0.5 * pairwise(SqEuclidean(), x_star'))
f_star = rand(MvNormal(K), 1)
fails with a Base.LinAlg.PosDefException. In this example, the points are close together and K becomes singular i.e. det(K) # 0.0.
The equivalent calculation in R:
x.star <- seq(-5, 5, by=0.1)
K = exp(-0.5 * rdist::pdist(x.star) ^ 2)
f = mvtnorm::rmvnorm(1, sigma=K)
works fine, because the mvtnorm package only requires Positive Semi-definite matrices for the covariance (due to the multiple options for matrix decomposition I guess? The “cholesky” method fails, but “svd” and “eigen” both work fine).
This leads me to my question: is it possible to sample from a multivariate normal distribution in Julia with a PSD matrix? Is there some argument/coercion that I’m missing?
By the way, you can “hack” your way around this error by adding a small constant to the diagonal elements of K, e.g.
x_star = -5:0.1:5
K_psd = exp(-0.5 * pairwise(SqEuclidean(), x_star'))
K_pd = copy(K_psd);
for i in 1:length(x_star)
K_pd[i, i] = K_pd[i, i] + 0.001
end
det(K_pd) # 2.00e-256
det(K_psd) # 0.00
f_star = rand(MvNormal(K_pd), 1)
but it would be nice to not have to do this.
