Transform random numbers into a valid covariance matrix

Question

I want to transform $n(n+1) / 2$ real random numbers into a valid covariance matrix $C \in \mathbb{R}^{n \times n}$. So far, I have tried to do this using the Cholesky decomposition of $C$, i.e. $C=L \cdot L^*$, where $L$ is a lower triangular matrix. So I basically populate the diagonal and the lower triangular part of $L$ with my random numbers and compute $C$. This works fine most of the time, however, I think due to numerical issues, I sometimes still end up with a matrix that is not positive-semidefinite.

Is there another, numerically safe way to achieve what I want to do? Or is there a good way to somehow regularize my "covariance" matrix in case it is not positive-semidefinite?

Edit (reply to Jan Kuckaka's answer)

If I understand correctly, you write the covariance matrix as $\Sigma = (\mathrm{diag}(\Sigma))^{1/2} \cdot \mathrm{corr} \cdot (\mathrm{diag}(\Sigma))^{1/2}$, where $\mathrm{corr}$ is the correlation matrix and $\mathrm{diag}(\Sigma)$ is a diagonal matrix holding the variances (see also Wikipedia article on covariance matrix). I tried this out using the following python code:

import numpy as np

n_dim = 3

std = np.random.uniform(low=0, high=5, size=n_dim)
corr = np.random.uniform(low=-1, high=1, size=int(n_dim * (n_dim - 1) / 2)) 

std_m = np.eye(n_dim) * std

corr_m = np.ones((n_dim, n_dim))
corr_m[np.triu_indices(n_dim, k=1)] = corr
corr_m[np.tril_indices(n_dim, k=-1)] = corr

cov_m = np.dot(np.dot(std_m, corr_m), std_m)

print(np.linalg.eigvals(cov_m))

I sample the standard deviations as random positive numbers and the correlations as random positive or negative numbers. For 2 dimensions, it seems to work, however, for higher dimensions, you can easily get negative eigenvalues. For example:

std = np.array([2.73023471, 4.83108501, 0.77834115])
corr = array([-0.65499313,  0.59416971,  0.61383891])

Edit Possible duplicate of this question

I think my question is not answered here. This question is about how one can sample random covariance matrices. My questions is: given $n(n+1)/2$ real numbers (wherever they might come from), is there a numerically safe way to form a valid covariance matrix from these numbers?

Turns out that 1) my answer is wrong. Thanks for pointing that out. 2) There is already an answer [here](https://stats.stackexchange.com/questions/2746/how-to-efficiently-generate-random-positive-semidefinite-correlation-matrices) — Jan Kukacka, Jun 21 '18 at 11:49
Could you please explain the sense in which you mean "transform"? Obviously it implies some rearrangement and duplication, but could it extend to other functions from $\mathbb{R}^{n(n+1)/2}$ into the the space of $n\times n$ positive semi-definite matrices? Which functions, exactly? — whuber, Jun 21 '18 at 12:05
Let me give some context to clarify: I am playing around with convolutional neural networks at the moment. I want to train a CNN to predict a covariance matrices. For $n$ dimensions, the network would predict $n(n+1)/2$ real numbers. Now I need to apply some transformation to these numbers in order to obtain a valid covariance matrix. So my requirement would be that the transformation is differentiable. The Cholesky approach almost does the job, however, numerical issues stop me from using it right now. — Joerg, Jun 21 '18 at 12:28
It sounds like you are looking for a *parameterization* of the manifold of covariance matrices. A natural one is to parameterize the orthogonal matrices separately from the singular values (and consider using the logarithms of the singular values for their parameters). But whether this is appropriate or not depends on your application and how you intend to measure predictive accuracy. — whuber, Jun 22 '18 at 13:55

score 0 · Answer 1 · answered Jun 21 '18 at 08:17

The notation that you use for your decomposition suggests that you compute C through the product of L and its conjugate transpose. Since you have only real numbers, you could just use the matrix transpose. Perhaps forcing the conjugate transpose leads to numerical issues with very small imaginary parts.

Transform random numbers into a valid covariance matrix

1 Answers1