# Transform random numbers into a valid covariance matrix [duplicate]

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?

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?

## marked as duplicate by Jan Kukacka, Ferdi, whuber♦Jun 22 '18 at 13:55

• 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