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

[thouis]

Original ticket http://projects.scipy.org/numpy/ticket/1223 Reported 2009-09-15 by trac user zero79, assigned to unknown.

numpy.random.multivariate_normal will happily accept non-physical covariance matrices (i.e. non-symetric and those with det(C)<0), which results in erroneous distributions. attached is code that demonstrates the problem for a covariance matrix with det(C)<0 (except when rhoxz=1, which produces a valid covariance matrix).

!/usr/bin/python

import numpy , pylab

nmcsamples = 210*6

xavg = 1.0 yavg = 2.0 zavg = 3.0

xdev = 0.1 ydev = 0.1 zdev = 0.1

rhoxy = 1.0 rhoxz = numpy.arange( 0.0 , 1.05 , 0.1 ) rhoyz = 1.0

nsamples = len( rhoxz ) xdevmc = numpy.zeros( nsamples ) ydevmc = numpy.zeros( nsamples ) zdevmc = numpy.zeros( nsamples ) rhoxymc = numpy.zeros( nsamples ) rhoxzmc = numpy.zeros( nsamples ) rhoyzmc = numpy.zeros( nsamples )

for n in range( 0 , nsamples ):

cxx = xdev**2
cyy = ydev**2
czz = zdev**2
cxy = rhoxy*xdev*ydev
cxz = rhoxz[n]*xdev*zdev
cyz = rhoyz*ydev*zdev

c = numpy.array( [ [ cxx , -cxy , cxz ] , [ cxy , cyy , cyz ] , [ cxz , cyz , czz ] ] )
x,y,z = numpy.transpose( numpy.random.multivariate_normal( [ xavg , yavg , zavg ] , c , nmcsamples ) )
print rhoxz[n],numpy.linalg.det( c )

xavgmc = numpy.mean( x )
yavgmc = numpy.mean( y )
zavgmc = numpy.mean( z )

xdevmc[n] = numpy.sqrt( numpy.mean( ( x - xavgmc )**2 ) )
ydevmc[n] = numpy.sqrt( numpy.mean( ( y - yavgmc )**2 ) )
zdevmc[n] = numpy.sqrt( numpy.mean( ( z - zavgmc )**2 ) )

rhoxymc[n]  = ( numpy.mean( x*y ) - xavgmc*yavgmc )/xdevmc[n]/ydevmc[n]
rhoxzmc[n]  = ( numpy.mean( x*z ) - xavgmc*zavgmc )/xdevmc[n]/zdevmc[n]
rhoyzmc[n]  = ( numpy.mean( y*z ) - yavgmc*zavgmc )/ydevmc[n]/zdevmc[n]

pylab.subplots_adjust( wspace = 0.4 ) pylab.rc( 'font' , size = 9 ) pylab.rc( 'legend' , fontsize = 9 )

pylab.subplot( 2 , 3 , 1 ) pylab.plot( rhoxz , numpy.ones( nsamples )xdev , rhoxz , xdevmc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\Delta x$' ) pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] ) pylab.subplot( 2 , 3 , 2 ) pylab.plot( rhoxz , numpy.ones( nsamples )ydev , rhoxz , ydevmc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\Delta y$' ) pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] ) pylab.subplot( 2 , 3 , 3 ) pylab.plot( rhoxz , numpy.ones( nsamples )zdev , rhoxz , zdevmc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\Delta z$' ) pylab.legend( ( 'input' , 'monte carlo' ) ) pylab.axis( [ 0.0 , 1.0 , 0.09 , 0.12 ] ) pylab.subplot( 2 , 3 , 4 ) pylab.plot( rhoxz , numpy.ones( nsamples )rhoxy , rhoxz , rhoxymc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\rho{xy}$' ) pylab.subplot( 2 , 3 , 5 ) pylab.plot( rhoxz , rhoxz , rhoxz , rhoxzmc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\rho{xz}$' ) pylab.subplot( 2 , 3 , 6 ) pylab.plot( rhoxz , numpy.ones( nsamples )*rhoyz , rhoxz , rhoyzmc ) pylab.xlabel( '$\rho{xz}$' ) pylab.ylabel( '$\rho_{yz}$' )

pylab.savefig( 'multivariate_normal-problem.pdf' )

[thouis]

Attachment in Trac by trac user zero79, 2009-09-15: multivariate_normal-problem

[thouis]

Comment in Trac by atmention:charris, 2009-09-15

''copy and paste formatting didn't work.''

You need to highlight the code and hit the code block button.

Your example code works fine here. If there is a negative singular value the svd code is likely the culprit and that indicates a problem with lapack/atlas. Where did you get those packages? What os/distro are you running? What numpy version?

[thouis]

Comment in Trac by trac user zero79, 2009-09-16

this code should execute just fine. the real problem is actually demonstrated in the plot file saved to multivariate_normal-problem.pdf.

this is not a problem with the svd code itself. its that non-physical covariances are permissable as inputs to multivariate_normal, where that should not be allowed. i've intentionally created bad covariance matrices (with negative determinants, and one with valid determinant, to demonstrate the problem). i would like to see numpy check for these cases and throw an error. otherwise, the user gets something that looks feasible, but is completely wrong.

fyi, this was tested on numpy 1.3.0 on debian lenny and sid.

[thouis]

Comment in Trac by trac user dgoldsmith, 2010-06-29

Let me be more succinct: the problem is not that the code doesn't run, the problem is that the code does run! Here's a much, much simpler example:

import numpy as N N.version.version '1.4.1' from numpy import random as R R.multivariate_normal((1,1), ((1,0),(0,-1))) # this should not run array([ 0.67651334, 0.44764747]) This shouldn't run because the covariance matrix is not "non-negative definite" (aka, positive semidefinite).

I discovered this bug independently while working on this function's docstring which states at the end of the Notes (i.e., almost as an afterthought):

"Note that the covariance matrix must be non-negative definite."

I thought it odd that such a requirement should be relegated to the Notes and not be included in the discussion of the cov function parameter, so I checked the code and found no "enforcement," so I checked the behavior, and found that, sure enough, the function is all too "happy" to accept input which the Notes state it "must" not be fed.

Now, I can see the reason - speed: my first call of this function w/ a two-component integer mean and a very simple integer-valued cov resulted in a brief but noticeable delay in being provided the answer; is this due to the use of svd, which a note in the code says should be replaced by Cholesky decomp.? - why we may want to leave it up to the user to enforce positive semidefiniteness (PSDness), but IMHO that's a dangerous road to go down, and if we do go down it, then the warning sign needs to be featured much more prominently, i.e., in the Extended summary and again in the cov Parameter description.

However, a better alternative, IMO, is to add code that checks for PSDness and raises a warning if it's detected - that should be the default behavior, but we can add a keyword parameter which allows the user to by-pass that check if they want to.

[thouis]

Comment in Trac by atmention:vincentdavis, 2010-06-29

I agree with this. (including adding the keyword) "However, a better alternative, IMO, is to add code that checks for PSDness and raises a warning if it's detected - that should be the default behavior, but we can add a keyword parameter which allows the user to by-pass that check if they want to. "

[thouis]

Comment in Trac by atmention:josef-pkt, 2010-06-29

some discussion is in the thread at

http://mail.scipy.org/pipermail/numpy-discussion/2009-September/045255.html

[thouis]

Comment in Trac by trac user dgoldsmith, 2010-06-30

Thanks, josef. Based on that discussion, might I suggest a hierarchical approach: first, Cholesky is tried to take the matrix square-root - according to CH in the above thread, IIUC, this will automatically provide the function w/ a check on the "physicality" of cov. If Cholesky works, great, we're done; if it throws an exception, to satisfy RK, an independent check on cov is performed (what would be the speediest way to independently assess cov for both symmetry and PSDness?). If cov passes that test, eigh (as suggested by CH) is tried; if it succeeds, the result is returned, along w/ a message that Cholesky failed, cov might be asymmetric and/or not PSD; if eigh also fails, give up: return np.empty() and a message that cov is likely asymmetric and/or not PSD.

[thouis]

Comment in Trac by atmention:mwiebe, 2011-03-24