Repository: scipy/scipy
Pull request: MAINT: stats: more accurate error message for `multivariate_normal`
#### Reference issue
Closes gh-15508

#### What does this implement/fix?
gh-15508 reported that `multivariate_normal` reported the error `singular matrix` with an input that was not singular. This PR makes the error message more accurate and adds notes to the documentation to explain that symmetry of the matrix is assumed. (gh-15579 explored checking for symmetry, but participants decided against adding the check.)


Problem statement:
BUG: Incorrect error message in multivariate_normal
### Describe your issue.

Hi everyone,

I'm getting an incorrect (?) error message when using multivariate normal. Here is my minimally reproducible example.

```python
from scipy.stats import multivariate_normal

cov = [[1., 0.],
       [1., 1.]]

f = multivariate_normal(cov=cov)
```
Output:
```bash
LinAlgError: singular matrix
```

This matrix has a non-zero determinant, is invertible, has rank 2, etc. It is clearly not a singular matrix, so I don't get where this error is coming from. 

As covariance matrices should (symmetric) positive definite (and this matrix is not), I guess the error in this case should be either:

* `not symmetric` or,
* `LinAlgError: Matrix is not positive definite`

Especially because one could use the `allow_singular=True` flag and `multivariate_normal` won't raise an error despite the cov matrix not being (symmetric) positive definite.

Thanks!

Ignacio

### Reproducing Code Example

```python
from scipy.stats import multivariate_normal

cov = [[1., 0.],
       [1., 1.]]

f = multivariate_normal(cov=cov)
```


### Error message

```shell
LinAlgError: singular matrix
```


### SciPy/NumPy/Python version information

1.6.2

Hints from issue discussion:
Hey @AtsushiSakai @IgnacioHeredia, I have submitted a PR to fix this error message. Please review it.
I'm not sure this is a bug or not, because a correct covariance matrix is always symmetric and positive semi-definite, so I think passing a non-symmetric matrix is a user error. 
@AtsushiSakai @NamamiShanker I glanced at gh-15579 but did not follow closely. If this error occurs for other non-symmetric matrices, why not just change the error message to "matrix is not symmetric positive definite" (or "semidefinite when `allow_singular` is true)? It seems like that is a simple change that would make the error message correct in strictly more cases than it is now. It is a user error, but it is still better if the error message were not inaccurate, right?
@mdhaber I’m sorry I don’t remember this discussion well… But i agree better error massage is always great.

Commit message: Update scipy/stats/_multivariate.py

[skip ci]

Co-authored-by: Tirth Patel <tirthasheshpatel@gmail.com>

Patch to review:
diff --git a/scipy/stats/_multivariate.py b/scipy/stats/_multivariate.py
index 025ffe2f924f..5fd1e674f8c1 100644
--- a/scipy/stats/_multivariate.py
+++ b/scipy/stats/_multivariate.py
@@ -158,10 +158,13 @@ def __init__(self, M, cond=None, rcond=None, lower=True,
 
         eps = _eigvalsh_to_eps(s, cond, rcond)
         if np.min(s) < -eps:
-            raise ValueError('the input matrix must be positive semidefinite')
+            msg = "The input matrix must be symmetric positive semidefinite."
+            raise ValueError(msg)
         d = s[s > eps]
         if len(d) < len(s) and not allow_singular:
-            raise np.linalg.LinAlgError('singular matrix')
+            msg = ("When `allow_singular is False`, the input matrix must be "
+                   "symmetric positive definite.")
+            raise np.linalg.LinAlgError(msg)
         s_pinv = _pinv_1d(s, eps)
         U = np.multiply(u, np.sqrt(s_pinv))
 
@@ -229,12 +232,12 @@ def random_state(self, seed):
 
 
 _mvn_doc_default_callparams = """\
-mean : array_like, optional
-    Mean of the distribution (default zero)
-cov : array_like, optional
-    Covariance matrix of the distribution (default one)
-allow_singular : bool, optional
-    Whether to allow a singular covariance matrix.  (Default: False)
+mean : array_like, default: ``[0]``
+    Mean of the distribution.
+cov : array_like, default: ``[1]``
+    Symmetric positive (semi)definite covariance matrix of the distribution.
+allow_singular : bool, default: ``False``
+    Whether to allow a singular covariance matrix.
 """
 
 _mvn_doc_callparams_note = """\
@@ -293,8 +296,11 @@ class multivariate_normal_gen(multi_rv_generic):
     -----
     %(_mvn_doc_callparams_note)s
 
-    The covariance matrix `cov` must be a (symmetric) positive
-    semi-definite matrix. The determinant and inverse of `cov` are computed
+    The covariance matrix `cov` must be a symmetric positive semidefinite
+    matrix when `allow_singular` is True; it must be (strictly) positive
+    definite when `allow_singular` is False.
+    Symmetry is not checked; only the lower triangular portion is used.
+    The determinant and inverse of `cov` are computed
     as the pseudo-determinant and pseudo-inverse, respectively, so
     that `cov` does not need to have full rank.
 
@@ -683,13 +689,13 @@ def __init__(self, mean=None, cov=1, allow_singular=False, seed=None,
 
         Parameters
         ----------
-        mean : array_like, optional
-            Mean of the distribution (default zero)
-        cov : array_like, optional
-            Covariance matrix of the distribution (default one)
-        allow_singular : bool, optional
-            If this flag is True then tolerate a singular
-            covariance matrix (default False).
+        mean : array_like, default: ``[0]``
+            Mean of the distribution.
+        cov : array_like, default: ``[1]``
+            Symmetric positive semidefinite covariance matrix of the
+            distribution.
+        allow_singular : bool, default: ``False``
+            Whether to allow a singular covariance matrix.
         seed : {None, int, `numpy.random.Generator`,
                 `numpy.random.RandomState`}, optional
 
diff --git a/scipy/stats/tests/test_multivariate.py b/scipy/stats/tests/test_multivariate.py
index 40841b268523..ce6ee7301c09 100644
--- a/scipy/stats/tests/test_multivariate.py
+++ b/scipy/stats/tests/test_multivariate.py
@@ -325,6 +325,13 @@ def test_exception_singular_cov(self):
         assert_raises(e, multivariate_normal.cdf, x, mean, cov)
         assert_raises(e, multivariate_normal.logcdf, x, mean, cov)
 
+        # Message used to be "singular matrix", but this is more accurate.
+        # See gh-15508
+        cov = [[1., 0.], [1., 1.]]
+        msg = "When `allow_singular is False`, the input matrix"
+        with pytest.raises(np.linalg.LinAlgError, match=msg):
+            multivariate_normal(cov=cov)
+
     def test_R_values(self):
         # Compare the multivariate pdf with some values precomputed
         # in R version 3.0.1 (2013-05-16) on Mac OS X 10.6.


Code review comments:
diff_hunk: @@ -683,13 +689,13 @@ def __init__(self, mean=None, cov=1, allow_singular=False, seed=None,
 
         Parameters
         ----------
-        mean : array_like, optional
-            Mean of the distribution (default zero)
-        cov : array_like, optional
-            Covariance matrix of the distribution (default one)
-        allow_singular : bool, optional
-            If this flag is True then tolerate a singular
-            covariance matrix (default False).
+        mean : array_like, default: ``[0]``
+            Mean of the distribution.
+        cov : array_like, default: ``[1]``
+            Symmetric positive semidefinite covariance matrix of the
path: scipy/stats/_multivariate.py
text: @user1:
```suggestion
            Symmetric positive (semi)definite covariance matrix of the
```

@user1:
Forgot to suggest for frozen distribution