[SciPy-user] How to generate random positive definite matrix

Robert Kern robert.kern@gmail....
Mon Jun 11 00:48:25 CDT 2007


David Cournapeau wrote:
> Hi there,
> 
>     I need to generate random positive definite matrix, mainly for 
> testing purpose. Before, I was generating them using a random matrix A 
> given by randn, and computing A'A. Unfortunately, if A is singular, so 
> is A'A. Is there a better way to do than testing whether A is singular ? 
> They do not need to follow a specific distribution (but I would like to 
> avoid them to follow a really special "pattern").

I would generate N random direction vectors (draw from a multivariate normal
distribution with eye(N) as the covariance matrix and normalize the samples to
be unit vectors). Resample any vector which happens to be nearly parallel to
another (i.e. the dot product is within some eps of 1). Now, form a correlation
matrix using the dot products of each of the unit vectors. Draw N random
positive values from some positive distribution like log-normal or gamma.
Multiply this vector on either side of the correlation matrix:

  v * corr * v[:,newaxis]

You now have a random positive definite matrix which is even somewhat interpretable.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco


More information about the SciPy-user mailing list