Constraints on generating valid covariance matrices for multivariate normal distributions

Question

I am interested in randomly generating multivariate normal distributions (MVND) as the underlying probability function to generate instances for a data stream. I understand that to do so requires two things (i.e. a MVND is parameterized by two things):

1. A vector in the $n$-dimensional space representing where the MVND is centered; and

2. An $n \times n$ matrix representing the covariance of the distribution.

Step 1 is trivial to generate and I have no problems there. Step 2, however, presents some challenges that I am trying to overcome. At the moment I am generating random values $\in [0.0, 1.0]$ to populate the covariance matrix, $C$. From what I have found online this process is subject to the following constraints:

1. $C_{ii}$ is equal to the variance of the $i$th random variable;

2. $C$ must be symmetric: $C_{ij} = C_{ji}$; and

3. $C$ must be positive definite: $z^TCz > 0 \text{ } \forall \text{ non-zero column vector }z\text{ of }n\text{ real numbers}$.

Constraint 3 seems to be the trickiest one to meet and was chosen so that the MVND's density function is "non-degenerate". I have been successful generating diagonal covariances matrices (described here) where $C_{ij} = C_{ji} = 0 \text{, }\forall i\neq j$, but this is a bit of a cheat. I would like to be able to fully simulate the full range of MVND and, eventually other multivariate distributions. I would like to know, then:

Are my constraints correct?

and

Am I missing any constraints?

Answer 1

Yes, your constraints are correct (though I wouldn't really call the first one a constraint, necessarily...).

One way to generate symmetric positive semi-definite matrices is as follows: if $X \in \mathbb R^{n \times p}$, then $X X^T \in \mathbb R^{n \times n}$ is psd: $(X X^T)^T = X X^T$ so that it's symmetric, and $z^T X X^T z = \lVert X^T z \rVert^2 \ge 0$. So any way that you choose $X$ will give you a valid covariance $X X^T$, and indeed any psd matrix can be decomposed in such a way (e.g. the Cholesky decomposition).

If you choose the entries of $X$ to be iid normal, for example, you get a Wishart distribution.

Answer 2

Your constraints are correct. This question was asked here. The most intuitive solution is sampling from an inverse Wishart distribution.

Answer 3

While this does not directly answer your question, the following provides an alternative procedure to produce random numbers according to a random MVND.

Consider $n$ independent random variables distributed according to the standard normal ($\mu=0$,$\sigma=1$), $G_j$. Your can form $n$ random variables $$X_i = \mu_i +\sum_{j=1}^n\,a_{ij}\,G_j$$ where $a_{ij}$ are constants. $X_i$ are distributed according to a MVND, with covariance, $$ V_{ij} = \sum_{k=1}^n a_{ik}\,a_{jk}$$ The particular MVND is defined by the choice of the elements of the matrix $A$. The full space of possible MVND is sampled with $A$ being an upper (or lower) triangle matrix.