BUG, API: np.random.multivariate_normal behavior with bad covariance matrix #5726

cowlicks · 2015-03-27T15:48:05Z

Edit:

This now allows for PSD checking with a check_valid keyword arg that can take the values 'warn', 'raise', and 'ignore'.

It also allows the covariance matrix to be almost PSD up to some small errors.

There have been several old discussions about how np.random.multivariate_normal should behave when passed a "bad" covariance matrix. i.e. non-symmetric, or not Positive semi-definite (PSD). The reason people need this is because the covariance matrix might have roundoff error making it not strictly symmetric or PSD.

In the discussions there was some consensus to fix things the way I do in this PR. I added a keyword arg for ignoring the PSD check and silencing the warning. I also added a keyword arg that ignores the new symmetry check and silences the warning. Previously symmetry was not checked at all.

I'd like to go one step further and raise ValueError when the covariance matrix isn't symmetric or PSD, of course, with the keyword args there so this can be bypassed. What do y'all think about that?

old discussions:
#1821
http://mail.scipy.org/pipermail/numpy-discussion/2009-September/045255.html

argriffing · 2015-03-27T16:19:38Z

Just as a technical comment, is there a reason to use svd instead of eigh here? The covariance matrix has already been constructed, and the line
neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)
looks weird to me. It seems like it should just directly check the approximate signs of eigenvalues instead, avoiding this matrix-matrix product.

cowlicks · 2015-03-27T16:28:45Z

The was a discussion for switching to eigh, or cholesky. eigh seems better. here but that seemed out of scope of this PR, but it is on my todo list.

In the code there is a comment that says We continue to use the SVD rather than Cholesky in order to preserve current outputs.. I'm guessing using eigh or Cholesky would change the results to some number of sig figs.

argriffing · 2015-03-27T16:36:36Z

OK, sorry for rehashing the old discussion.

charris · 2015-03-27T17:20:15Z

Needs to be discussed on the list. Another option is to add a method keyword and an eigh method that checks for symmetry by default. The semidefinite property is not well defined for floating point for small eigenvalues due to roundoff, so there is the question of what to do in that case. It would also probably be easier to add a single check keyword with boolean values and either go as is or check both symmetry and positive semidefinite.

cowlicks · 2015-03-27T19:34:12Z

@charris I'll ping the list.

On Travis-CI a RuntimeWarning isn't being caught even though it is gaurded by with warings.catch_warnings .... Failure here, and guarding code here.

This doesn't happen on my local machine, what is different on Travis-CI?

charris · 2015-03-27T20:12:48Z

Try warnings.simplefilter("always") as the first line in the managed context.

charris · 2015-03-27T20:24:08Z

Note that symettry and almost positive semidefinite can be checked with the svd. Something like

In [16]: a = eye(3)

In [17]: u, d, v = svd(a)

In [18]: b = dot(v.T*d, v)

In [19]: allclose(a, b)
Out[19]: True

The tolerances for allclose can be adjusted.

cowlicks · 2015-04-10T15:24:31Z

@charris It didn't generate much discussion on list (link), but I think @josef-pkt's comments were helpful. To summarize, instead of seperate kwargs for checking symmetry, and PSDness, just one kwarg check_valid= that takes 'warn', 'raise', and 'ignore'. And it will only assert symmetry up to some reasonable tolerance (instead of exact).

cowlicks · 2015-05-10T21:44:46Z

@charris Rehashed this. Basically we now allow for small floating point errors in the covariance matrix. And I added a check_valid key word that lets you choose to raise, warn, or ignore if the covariance matrix is PSD.

cowlicks · 2015-05-10T21:45:22Z

ping @josef-pkt

josef-pkt · 2015-05-10T23:28:33Z

numpy/random/mtrand/mtrand.pyx

        (u, s, v) = svd(cov)
-        neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)
+


I would be the entire next part under
if check_valid != 'ignore':

no reason to check if we 'ignore'

In np.random.multivariate_normal allow the covariance matrix to have small floating point errors. And allow control over what to do if the PSD check fails.

Test that warnings are or are not raised when check_valid kwarg is used.

homu · 2016-01-24T00:20:18Z

☔ The latest upstream changes (presumably #7082) made this pull request unmergeable. Please resolve the merge conflicts.

rgommers · 2017-01-03T20:35:08Z

@ovillellas pushed 0 commits.

You managed to confuse GitHub, haven't seen that before:)

Have you taken over this PR, and is it ready to merge?

ovillellas · 2017-01-03T20:36:48Z

Yes, I just took over this PR. It should be ready to commit (it was just a matter of cleaning up a bit the test).

rgommers · 2017-01-03T20:39:40Z

Okay thanks. I'll review it now then, before we forget about it again.

rgommers · 2017-01-03T20:44:38Z

numpy/random/mtrand/mtrand.pyx

@@ -4380,6 +4381,10 @@ cdef class RandomState:
            generated, and packed in an `m`-by-`n`-by-`k` arrangement.  Because
            each sample is `N`-dimensional, the output shape is ``(m,n,k,N)``.
            If no shape is specified, a single (`N`-D) sample is returned.
+        check_valid : 'warn', 'raise', 'ignore'


formatting: change to {'warn', 'raise', 'ignore'}, optional

rgommers · 2017-01-03T20:45:07Z

numpy/random/mtrand/mtrand.pyx

@@ -4380,6 +4381,10 @@ cdef class RandomState:
            generated, and packed in an `m`-by-`n`-by-`k` arrangement.  Because
            each sample is `N`-dimensional, the output shape is ``(m,n,k,N)``.
            If no shape is specified, a single (`N`-D) sample is returned.
+        check_valid : 'warn', 'raise', 'ignore'
+            Behavior when the covariance matrix is not Positive Semi-definite.


style nitpick: use lower case for positive semidefinite

rgommers · 2017-01-03T20:46:05Z

numpy/random/mtrand/mtrand.pyx

@@ -4380,6 +4381,10 @@ cdef class RandomState:
            generated, and packed in an `m`-by-`n`-by-`k` arrangement.  Because
            each sample is `N`-dimensional, the output shape is ``(m,n,k,N)``.
            If no shape is specified, a single (`N`-D) sample is returned.
+        check_valid : 'warn', 'raise', 'ignore'
+            Behavior when the covariance matrix is not Positive Semi-definite.
+        tol : float


float, optional

rgommers · 2017-01-03T20:49:51Z

numpy/random/mtrand/mtrand.pyx

        if mean.shape[0] != cov.shape[0]:
-               raise ValueError("mean and cov must have same length")
+            raise ValueError("mean and cov must have same length")



Can you add a check for check_valid being one of the 3 accepted strings?

actually, move it from line 4511 up here. It's misplaced there; only checked for when not psd.

that looks good now

note that the check is itself later in the code, when checking. However, it may go unchecked when the psd passes, so I have rearranged the logic so it is always checked.

rgommers · 2017-01-03T20:52:56Z

numpy/random/mtrand/mtrand.pyx

+        check_valid : 'warn', 'raise', 'ignore'
+            Behavior when the covariance matrix is not Positive Semi-definite.
+        tol : float
+            Tolerance of the singular values in covariance matrix.


absolute or relative tolerance?

rgommers · 2017-01-03T20:53:20Z

numpy/random/mtrand/mtrand.pyx

-                          RuntimeWarning)
+
+        if check_valid != 'ignore':
+            psd = np.allclose(np.dot(v.T * s, v), cov)


tol should be used here. It's not used anywhere now.

Can you add a test with non-default tol for good measure?

I realize. I wonder if tol should be changed to atol and rtol and just forward them in here.

You'll have to specify both atol and rtol (here, not in the signature), because the defaults are both nonzero: rtol=1e-05, atol=1e-08

Sorry my comment wasn't a reply to yours, non-updating interface.

atol = rtol = tol seems okay. atol is the only one of interest near covariance values of 0, and rtol may do something useful for very large values.

rgommers · 2017-01-03T21:04:32Z

Overall LGTM, let's get this merged once comments are addressed.

ovillellas · 2017-01-03T21:07:42Z

Ok, just pushed the documentation changes.

Also:

changed logic handling check_valid so it will always raise on invalid argument even if it cov is not psd.
changed tol parameter to atol and rtol, that are mapped straight into np.allclose in the code.

updated the docstring accordingly

10000

rgommers · 2017-01-03T21:13:19Z

changed tol parameter to atol and rtol, that are mapped straight into np.allclose in the code.

I wonder if that's better. Realistically no one will use rtol for anything interesting, so likely adding only one keyword is better. What do you think about atol = rtol = tol?

ovillellas · 2017-01-03T21:19:48Z

Ok, I will change it back to the single tol and pass it as atol and rtol

charris · 2017-01-03T21:25:34Z

Needs a mention in doc/release/1.13.0-notes.rst.

ovillellas · 2017-01-03T21:37:57Z

reverted to a single tol argument. I also added a mention in the release notes under "changes"

rgommers · 2017-01-04T03:28:59Z

Okay, in it goes! Thanks @ovillellas, @cowlicks, all.

cowlicks force-pushed the mult-norm branch from e6d6ec3 to 86b729d Compare March 27, 2015 15:50

charris added 01 - Enhancement component: numpy.random labels Mar 28, 2015

cowlicks force-pushed the mult-norm branch from 57a43c6 to da8ae9f Compare May 10, 2015 21:42

josef-pkt reviewed May 10, 2015
View reviewed changes

cowlicks added 2 commits May 19, 2015 16:21

BUG, API: Allow covariance matrix with small fp errors.

fc57915

In np.random.multivariate_normal allow the covariance matrix to have small floating point errors. And allow control over what to do if the PSD check fails.

TST: test multivariate_normal check_valid kw

f555826

Test that warnings are or are not raised when check_valid kwarg is used.

cowlicks force-pushed the mult-norm branch from da8ae9f to f555826 Compare May 19, 2015 21:29

Oscar Villellas added 2 commits January 3, 2017 20:28

Merge remote-tracking branch 'numpy-org/master' into mult-norm

4c93e28

fixed merged test

fde2617

rgommers reviewed Jan 3, 2017

View reviewed changes

Documentation fix and proper handling of tolerance

6d7f14f

single too argument + mention in release docs.

c85d199

rgommers merged commit b94c2b0 into numpy:master Jan 4, 2017

rgommers added this to the 1.13.0 release milestone Jan 4, 2017

rgommers mentioned this pull request Jan 4, 2017

numpy.random.multivariate_normal accepts indefinite covariance matrices (Trac #1223) #1821

Closed

homu mentioned this pull request Jan 4, 2017

ENH: random: Add multivariate_hypergeometric function. #8056

Closed

rgommers mentioned this pull request Jan 5, 2017

It would be great to have multivariate_normal raising exception #8263

Closed

xieliaing mentioned this pull request Nov 7, 2017

MultivariateNormal can't handle semi-positive definite covariance matrix accord-net/framework#1019

Open

3 tasks

		(u, s, v) = svd(cov)
		neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)

Uh oh!

BUG, API: np.random.multivariate_normal behavior with bad covariance matrix #5726

BUG, API: np.random.multivariate_normal behavior with bad covariance matrix #5726

Uh oh!

Conversation

Edit:

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!