Description
instead of undefined behaviour as the NumPy doc says:
Note that the covariance matrix must be positive semidefinite (a.k.a. nonnegative-definite). Otherwise, the behavior of this method is undefined and backwards compatibility is not guaranteed.
For example, the following
from numpy import errstate
from numpy.random import RandomState
from numpy.linalg import eigh
from numpy import dot
rs = RandomState(3)
m = rs.randn(10)
G = rs.randn(10, 5)
K = dot(G, G.T)
with errstate(all='raise'):
try:
rs.multivariate_normal(m, K)
except Exception as e:
print("Exception raised!")
else:
print("Exception not raised!")
prints
bug.py:13: RuntimeWarning: covariance is not positive-semidefinite.
rs.multivariate_normal(m, K)
Exception not raised!
As you know, K
happens to not be positive-semidefinite in the above example because of floating point arithmetic. This can be easily fixed by adding a small value to its diagonal. If NumPy raised an Exception, I could perform this hack only if necessary. As it is now, I would need to check for positive-semidefiniteness beforehand, or add a small value to its diagonal before hand. Both solutions are not ideal in my opinion. Raising an exception would solve that.