Scipy.stats.multivariate_normal raising `LinAlgError: special matrix`, although my covariance matrix is ​​reversible

Question

I'm having trouble trying to use scipy.stats.multivariate_normal , hopefully one of you can help.

I have a 2x2 matrix that can find the reverse use of numpy.linalg.inv() , however, when I try to use it as the covariance matrix in multivariate_normal , I get a LinAlgError , stating that it is a special matrix:

 In [89]: cov = np.array([[3.2e5**2, 3.2e5*0.103*-0.459],[3.2e5*0.103*-0.459, 0.103**2]]) In [90]: np.linalg.inv(cov) Out[90]: array([[ 1.23722158e-11, 1.76430200e-05], [ 1.76430200e-05, 1.19418880e+02]]) In [91]: multivariate_normal([0,0], cov) --------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) <ipython-input-91-44a6625beda5> in <module>() ----> 1 multivariate_normal([0,0], cov) /mnt/ssd/Enthought_jli199/Canopy_64bit/User/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc in __call__(self, mean, cov, allow_singular, seed) 421 return multivariate_normal_frozen(mean, cov, 422 allow_singular=allow_singular, --> 423 seed=seed) 424 425 def _logpdf(self, x, mean, prec_U, log_det_cov, rank): /mnt/ssd/Enthought_jli199/Canopy_64bit/User/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc in __init__(self, mean, cov, allow_singular, seed) 591 """ 592 self.dim, self.mean, self.cov = _process_parameters(None, mean, cov) --> 593 self.cov_info = _PSD(self.cov, allow_singular=allow_singular) 594 self._dist = multivariate_normal_gen(seed) 595 /mnt/ssd/Enthought_jli199/Canopy_64bit/User/lib/python2.7/site-packages/scipy/stats/_multivariate.pyc in __init__(self, M, cond, rcond, lower, check_finite, allow_singular) 217 d = s[s > eps] 218 if len(d) < len(s) and not allow_singular: --> 219 raise np.linalg.LinAlgError('singular matrix') 220 s_pinv = _pinv_1d(s, eps) 221 U = np.multiply(u, np.sqrt(s_pinv)) LinAlgError: singular matrix 

Answer 1

By default, multivariate_normal checks if any of the eigenvalues ​​of the covariance matrix is ​​smaller than a specific tolerance selected based on its dtype type and the value of its largest eigenvalue (see source code for scipy.stats._multivariate._PSD and scipy.stats._multivariate._eigvalsh_to_eps for full details).

As @kazemakase mentioned above, while your covariance matrix may be reversible according to the criteria used by np.linalg.inv , it is still very poorly prepared and does not perform the more stringent tests used by multivariate_normal .

You can go through allow_singular=True to multivariate_normal to skip this test, but in general it would be better to rescale your data to avoid passing such a poorly conditioned covariance matrix.