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.