Numerically unstable Cholesky?


I’m trying to sample from a simple Gaussian process (without a diagonal jitter term). It seems that torch’s and numpy’s multivariate Gaussian implementations are slightly different:

import numpy as np, torch
dx = 0.1
N = 10
x = torch.arange(N).view([-1,1]) * dx
mean = torch.zeros(N)
dist = torch.cdist(x, x)
K = torch.exp(-0.5 * dist**2)
s = torch.distributions.MultivariateNormal(mean,K).sample() # error: cholesky_cpu: U(8,8) is zero, singular U.
s = np.random.multivariate_normal(mean.numpy(), K.numpy()) # executes

Both numpy and torch arrays are float32.

I’m using anaconda3 2019.07 on my Mac 10.12.6.
Python 3.7.3
Numpy 1.16.4
Torch 1.4.0