Is there some good multivariate gaussian generator?

1

Dear all, a quick question,

Is there some good multivariate gaussian generator?
Or, like, what multivariate gaussian generator do you use?
You know, it should be able to generate gaussian random vector with covariance or something.

Key word
multivariate

Is Alan Miller’s random.f90 OK for multivariate gaussian generator?
https://jblevins.org/mirror/amiller/random.f90

Alan Miller and John Burkardt are two big names! Like Tom and Jerry!

Thank you very much in advance!

2

I think this algorithm works directly with the covariance matrix (as input argument) which makes it nice and rather easy to work with, but not the most efficient for repeated random-number generations from the same distribution. As I said before in another comment, the most efficient method is to pass the Cholesky factorization to the procedure rather than passing the covariance matrix (which would require the factorization computation on each and every call to the procedure). Here is an implementation that takes the Cholesky Factorization instead of the covariance matrix.

3 Likes
3

You are probably going to be asking many such similar questions. I am going to put in a link to a NAG Library routine here

https://www.nag.com/numeric/nl/nagdoc_latest/flhtml/g05/g05rzf.html

One advantage of using a tested numerical library is that your program will not be an aggregate of codes of variable quality, provenance and robustness. The source you (or others) linked to could be impeccable but where is the test program that gives you confidence that the compiler, with the custom compiler options that you throw at it, did not mess it up? Our makefiles at NAG are littered with special cases for which reduced optimization is necessary for some versions of each compiler. I think readers of this thread should be aware of such issues.

3 Likes
4

Thank you very much @shahmoradi , I always remember this thread.
A quick naive question, does the covariance matrix has to be positive definite for Cholesky factorization?
If the covariance matrix is not positive definite, will that cause some problem? If so is there some solutions? Thanks much!

5

You’re welcome. Yes. A probability distribution’s covariance matrix is always positive-definite. The Cholesky factorization is, in fact, the best way of testing for the positive-definiteness of a matrix as far as I know. The Covariance matrix is an entity that determines the scale of the distribution along each dimension. If it is not positive-definite, then it has fewer dimensions than the matrix or data indicates. The model is likely wrong, or there are deterministic internal dependencies in the data.

1 Like
6

Awesome! Thank you! :100:

7

BTW: Covariance matrices are only guaranteed to be (symmetric) positive semi-definite as far as I’m aware…

2 Likes
8

Yeah. I think what @shahmoradi means is that, a meaningful covariance matrix should be positive definite. If the covariance matrix is not positive definite, it tells us some extra information, like, some variables may be the linear combination of some variables,

normal distribution - What does a non positive definite covariance matrix tell me about my data? - Cross Validated.

that makes the covariance matrix not full rank and therefore not positive-definite, that will cause Cholesky factorization error. if that happens, need to reduce the number of variables (therefore the corresponding model needs to be changed as well), then the new covariance matrix should be positive definite.

1 Like
9

For the case when the correlation matrix computed in the usual way is not positive definite, NAG has a routine that calculates the nearest correlation matrix. The paper The Nearest Correlation Matrix Problem: Solution by Differential Evolution Method of Global Optimization says that Fortran code by the author, Sudhanshu K. Mishra, is available on request.

2 Likes
10

Thanks @Beliavsky.
Just curious, is there a subroutine that can do an approximated Cholesky factorization when an non positive definite matrix is encountered?
Like

A=chol(H)

A is what we want, a lower triangle matrix.
However if H happen to be non positive definite matrix, is there things like ‘chol_approximate’,

A_approximate = chol_approximate(H)

which can find an approximated lower triangle matrix A_approximate?

But anyway, of course, this is just a workaround. The real fix is to make the H Hermitian and positive-definite. Thanks!

EDIT.
A workaorund can be found below,

Use eigenvalue decomposition to get the R matrix which approximate the L=Chol(H) matrix.

image
image1403×1477 307 KB