Confusing behavior of np.random.multivariate_normal

Question

I am sampling from a multivariate normal using numpy as follows.

mu = [0, 0]
cov = np.array([[1, 0.5], [0.5, 1]]).astype(np.float32)
np.random.multivariate_normal(mu, cov)

It gives me the following warning.

RuntimeWarning: covariance is not positive-semidefinite.

The matrix is clearly PSD. However, when I use a np.float64 array, it works fine. I need the covariance matrix to be np.float32. What am I doing wrong?

Answer 1

This has been fixed in March 2019. If you still see the warning consider updating your numpy.

The warning is raised even for very small off-diagonal elements > 0. The default tolerance value does not seem to work well for 32 bit floats.

As a workaround pass a higher tolerance to the function:

np.random.multivariate_normal(mu, cov, tol=1e-6)

---

Details

np.random.multivariate_normal checks if the covariance is PSD by first decomposing it with (u, s, v) = svd(cov), and then checking if the reconstruction np.dot(v.T * s, v) is close enough to the original cov.

With float32 the result of the reconstruction is further off than the default tolerance of 1e-8 allows, and the function raises a warning.