[Numpy-discussion] Issue with numpy.random.multivariate_normal Linux RHEL4 np version: 1.6.1

John Gleeson jdgleeson at mac.com
Sat May 12 03:05:00 EDT 2012

On 2012-05-11, at 4:01 PM, Norman Shelley wrote:

> Running on Linux RHEL4
> numpy.random.multivariate_normal seems to work well for 25 mean  
> values but when I go to 26 mean values (and their corresponding  
> covariance values) it produces garbage.
> Any ideas?

The implementation of multivariate_normal in numpy uses Singular Value  
Decomposition, and your
covariance matrices are very ill-conditioned. That is true, however,  
for both the good case and the bad case, so it really doesn't explain  
why the SVD method suddenly fails when the dimension changes from 25  
to 26. I hope to look into this when I get more time.

I created a new version of multivariate_normal (see below) that has a  
method parameter. Setting method='chol' chooses Cholesky decomposition  
instead of SVD. It gives nice results for your size 26 case.

data = multivariate_normal(means, covariances, 2, method='chol')
print 'using chol:'
for i, row in enumerate(data):
     print i, row

def multivariate_normal(mean, cov, size=None, method='svd'):
     mean = np.array(mean)
     cov = np.array(cov)
     if size is None:
         shape = []
         shape = size
     if len(mean.shape) != 1:
            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")
     if mean.shape[0] != cov.shape[0]:
            raise ValueError("mean and cov must have same length")
     # Compute shape of output
     if isinstance(shape, int):
         shape = [shape]
     final_shape = list(shape[:])
     # Create a matrix of independent standard normally distributed  
     # numbers. The matrix has rows with the same length as mean and as
     # many rows are necessary to form a matrix of shape final_shape.
     x = np.random.standard_normal(np.multiply.reduce(final_shape))
     x.shape = (np.multiply.reduce(final_shape[0:len(final_shape)-1]),
     # Transform matrix of standard normals into matrix where each row
     # contains multivariate normals with the desired covariance.
     # Compute A such that dot(transpose(A),A) == cov.
     # Then the matrix products of the rows of x and A has the desired
     # covariance.

     if method == 'svd':
         (u,s,v) = np.linalg.svd(cov)
         x = np.dot(x*np.sqrt(s),v)
     elif method == 'chol':
         L = np.linalg.cholesky(cov)
         x = np.dot(x,L)
         raise ValueError("method must be 'svd' or 'chol'")
     # The rows of x now have the correct covariance but mean 0. Add
     # mean to each row. Then each row will have mean mean.
     x.shape = tuple(final_shape)
     return x

More information about the NumPy-Discussion mailing list