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
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.