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

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.

1 files changed, 21 insertions, 9 deletions

diff --git a/numpy/random/mtrand/mtrand.pyx b/numpy/random/mtrand/mtrand.pyx
index 5c949d07b..32f268ad3 100644
--- a/numpy/random/mtrand/mtrand.pyx
+++ b/numpy/random/mtrand/mtrand.pyx
@@ -4240,7 +4240,8 @@ cdef class RandomState:
                            self.lock)
 
     # Multivariate distributions:
-    def multivariate_normal(self, mean, cov, size=None):
+    def multivariate_normal(self, mean, cov, size=None, check_valid='warn',
+                            tol=1e-8):
         """
         multivariate_normal(mean, cov[, size])
 
@@ -4265,6 +4266,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
+            Tolerance of the singular values in covariance matrix.
 
         Returns
         -------
@@ -4349,11 +4354,11 @@ cdef class RandomState:
             shape = size
 
         if len(mean.shape) != 1:
-               raise ValueError("mean must be 1 dimensional")
+            raise ValueError("mean must be 1 dimensional")
         if (len(cov.shape) != 2) or (cov.shape[0] != cov.shape[1]):
-               raise ValueError("cov must be 2 dimensional and square")
+            raise ValueError("cov must be 2 dimensional and square")
         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")
 
         # Compute shape of output and create a matrix of independent
         # standard normally distributed random numbers. The matrix has rows
@@ -4376,12 +4381,19 @@ cdef class RandomState:
         # not zero. We continue to use the SVD rather than Cholesky in
         # order to preserve current outputs. Note that symmetry has not
         # been checked.
+
         (u, s, v) = svd(cov)
-        neg = (np.sum(u.T * v, axis=1) < 0) & (s > 0)
-        if np.any(neg):
-            s[neg] = 0.
-            warnings.warn("covariance is not positive-semidefinite.",
-                          RuntimeWarning)
+
+        if check_valid != 'ignore':
+            psd = np.allclose(np.dot(v.T * s, v), cov)
+            if not psd:
+                if check_valid == 'warn':
+                    warnings.warn("covariance is not positive-semidefinite.",
+                                  RuntimeWarning)
+                elif check_valid == 'raise':
+                    raise ValueError("covariance is not positive-semidefinite.")
+                else:
+                    raise ValueError("check_valid must equal 'warn', 'raise', or 'ignore'")
 
         x = np.dot(x, np.sqrt(s)[:, None] * v)
         x += mean