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

[David Cournapeau]

Anne Archibald wrote:
> On 10/06/07, David Cournapeau <david at ar.media.kyoto-u.ac.jp> 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").
>
> Every symmetric positive definite matrix (you didn't say symmetric but
> I assume you want it...)
Yes, I need them to as parameters of multivariate gaussian. so they need 
to be real, not hermitian.
>  can be diagonalized using an orthogonal
> matrix. That is, you can always factor it as O'DO, where D is diagonal
> and has positive elements on the diagonal. So to get a random one, do
> it backwards: pick a set of eigenvalues from some distribution, then
> for each pair of axes, rotate the matrix by a random angle in the
> plane they form (which you can do with a cos/sin/-sin/cos matrix) with
> a randomly-decided reflection. This gives a distribution that's
> invariant under orthonormal change of basis; your only nonuniformities
> come in the choice of eigenvalues, and those you will probably want to
> control in testing. Eigenvalues of wildly differing sizes tend to lead
> to big numerical headaches, so you can experiment with the size of
> headache you want to give your code. If you want them all roughly the
> same size you could take normal variates and square them, but for
> better control I'd go with exponentiating a normal distribution.
That's a pretty good idea, actually. I totally forgot about the fact 
that the diagonalization vectors of an hermitian matrix can be 
interpreted as rotation (I really have to check my knowledge on those 
special matrices, I got all rusty, and I studied linear algebra less 
than 5 years ago :) ).
>
> There's an alternative, probably faster, approach based on the
> uniqueness of the Cholesky square root of a matrix - just choose your
> eigenvalues, again, put their square roots on the diagonal, and
> randomly generate the upper triangle. Then multiply A'A and you should
> be able to get an arbitrary (symmetric) positive definite matrix.
> Unless you are cleverer than I am about choosing the distributions of
> the entries, it won't be nicely direction-independent.
That's why I tried, but I found they are really "too special".
>
> Now, a sort-of interesting question is, is there a natural
> distribution on the cone of positive definite matrices one could hope
> to draw from? 
Well, we could draw samples from  a Wishart, but I was too lazy to 
implement it. It is definitely "natural", specially for multivariate 
normal, but I need to improve my knowledge on multidimensional calculus 
before being able to understand them (and implement them, hopefully).

cheers,

David