Invalid np.random.multivariate_normal behavior

Question

I take a sample from a multidimensional 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) 

He gives me the following warning.

RuntimeWarning: covariance is not positive-semidefinite.

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

Answer 1

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

As a workaround, skip the higher function tolerance:

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

---

More 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(vT * s, v) is close enough to the original cov .

With float32 the reconstruction result is even larger than the default allowable allowable 1e-8 , and the function raises a warning.