[Numpy-discussion] Behavior of np.random.multivariate_normal with bad covariance matrices

josef.pktd at gmail.com josef.pktd at gmail.com
Mon Mar 30 09:34:38 EDT 2015

On Sun, Mar 29, 2015 at 7:39 PM, Blake Griffith
<blake.a.griffith at gmail.com> wrote:
> I have an open PR which lets users control the checks on the input
> covariance matrix. The matrix is required to be symmetric and positve
> semi-definite (PSD). The current behavior is that NumPy raises a warning if
> the matrix is not PSD, and does not even check for symmetry.
> I added a symmetry check, which raises a warning when the input is not
> symmetric. And added two keyword args which users can use to turn off the
> checks/warnings when the matrix is ill formed. So this would only cause
> another new warning to be raised in existing code.
> This is needed because sometimes the covariance matrix is only *almost*
> symmetric or PSD due to roundoff error.
> Thoughts?

My only question is why is **exact** symmetry relevant?

A empirical covariance matrix might not be exactly symmetric unless we
specifically force it to be. But I don't see why some roundoff errors
that violate symmetry should be relevant.

use allclose with floating point rtol or equivalent?

Some user code might suddenly get irrelevant warnings.

neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)
doesn't need to be calculated if cov_psd is false.


some more:

svd can hang if the values are not finite, i.e. nan or infs

counter proposal would be to add a `check_valid` keyword with option
ignore. warn, raise, and "fix"

and raise an error if there are nans and check_valid is not ignore.


np.random.multivariate_normal   is only relevant if you have a new cov
each call (or don't mind repeated possibly expensive calculations),
so, I guess, adding checks by default won't upset many users.


> PR: https://github.com/numpy/numpy/pull/5726
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

More information about the NumPy-Discussion mailing list