I am trying to create random covariance matrices for three joint gaussian variables. My goal is to sample random covariances matrices that always have correlation between 0.7 and 0.9 (or 0 if there isn’t).

So far I am doing it manually with a repeat until is.positive.definite is true… But I am unable to achieve it, my repeat takes a lot of time because most of my matrices samples return false for the positive.definite.

Is there a library to do this or an simpler approach for this?

On the math side I know I can have correlation between: $X_1$ and $X_2$. $X_2$ and $X_3$. $X_1$ and $X_3$ If I am not mistaken, I can have correlation between the three pair or just one pair, there shouldn’t be any issue. But if there is correlation between two of them, the remaining correlation couldn’t be 0, otherwise the matrix would never be positive definite…

  • 5
    $\begingroup$ Covariance matrices just like vectors can be random variables with arbitrary distributions, so you cannot generate a "random" matrix without first specifying its distribution.The most common distribution is the Wishart distribution. $\endgroup$
    – Xi'an
    Commented Apr 12, 2020 at 3:58
  • 2
    $\begingroup$ See stats.stackexchange.com/questions/215497/…. There are many ways to adapt the ideas there. For instance, you could generate a sample from a trivariate Normal distribution having the desired underlying correlations: for sufficiently large samples, the sample covariance matrix will be very likely to satisfy your conditions. $\endgroup$
    – whuber
    Commented Apr 12, 2020 at 13:43

4 Answers 4


In this particular case there's a simple, easy, completely general method.

Break the problem down into two parts:

  1. Generate random variances $\sigma_i^2,$ $i=1,2,3.$ These define a diagonal matrix $\Sigma = \pmatrix{\sigma_1&0&0\\0&\sigma_2&0\\0&0&\sigma_3}.$

  2. Generate a random correlation matrix $R = \pmatrix{1&\rho_3&\rho_2\\\rho_3&1&\rho_1\\\rho_2&\rho_1&1}.$

The resulting random covariance is $\Sigma R \Sigma.$ It is symmetric by construction. It will be positive-definite if and only if $R$ is, which is equivalent to $|\rho_3|\le 1,$ $|\rho_2|\le 1,$ and $R$ has positive determinant.

What happens if you generate $(\rho_1,\rho_2,\rho_3)$ using any distribution you like supported on the cube $[0.7,0.9]^3$? The only condition you need to check concerns the determinant. But since

$$\det R = 1 - (\rho_1^2+\rho_2^2+\rho_3^2) + 2\rho_1\rho_2\rho_3,$$

we may do a little bit of Calculus and establish that the minimum value of the determinant is attained when one of the $\rho_i$ equals $0.7$ and the other two equal $0.9,$ with a value of $24/1000\gt 0.$ Consequently

no matter how $\rho_1, \rho_2, \rho_3$ are generated, $\det R$ is always positive. Therefore, provided the $\sigma_i$ are positive, $\Sigma R \Sigma$ is a positive-definite covariance matrix.

As an example, you could generate the $\sigma_i^2$ independently with (say) some Gamma distribution and generate the $\rho_i$ uniformly. I created $100,000$ such covariance matrices this way; it took less than two seconds. Here's a summary of the results on which are superimposed the intended distribution density functions, showing the method works as intended.


It is clear that

When $\sigma_1, \ldots, \rho_3$ are drawn from any six-dimensional distribution supported on $(0,\infty)^3\times (0.7,0.9)^3,$ $\Sigma R \Sigma$ is a valid covariance matrix with all correlations between $0.7$ and $0.9.$ Conversely, any distribution of covariance matrices with these properties determines such a distribution of $\sigma_1, \ldots, \rho_3.$

You can even introduce dependencies between the $\sigma_i$ and the $\rho_j$ if you like.

This is the R code to reproduce the figure. rcov generates an array of n such covariance matrices (referenced by a third index).

rcov <- function(n=1, shape=1, rate=1) {
  sigma <- matrix(rgamma(3*n, shape, rate), 3)
  rho <- matrix(runif(3*n, 0.7, 0.9), 3)
  array(sapply(1:n, function(i) {
    diag(sigma[,i]) %*% matrix(c(1, rho[3,i], rho[2,i],
                                rho[3,i], 1, rho[1,i],
                                rho[2,i], rho[1,i], 1), 3, 3) %*% diag(sigma[,i])
  }), c(3,3,n))

shape <- c(2, 5, 10)
rate <- shape
system.time(rho <- apply(Sigma <- rcov(1e5, shape, rate), 3, cov2cor)[c(2, 3, 6), ])

gray <- "#f0f0f0"
hist(rho, freq=FALSE, col=gray,
     main=expression(paste("Histogram of all ", rho[i])), xlab="Value")
abline(h=1 / (0.9 - 0.7), lwd=2)
for (i in 1:3) {
  hist(sqrt(Sigma[i,i,]), freq=FALSE, breaks=30, col=gray,
       main=bquote(sigma[.(i)]), xlab="Value")
  curve(dgamma(x, shape[i], rate[i]), lwd=2, add=TRUE)
  • $\begingroup$ Please fix the typo for $|\rho_3|\leq1$ and $|\rho_2|\leq1$ (it should be $<$, Sylvester's criterion). $\endgroup$
    – user289381
    Commented Jul 15, 2020 at 22:21
  • $\begingroup$ @ping That's no typo: the strict positivity of $\det(R)$ rules out the exceptions you are worried about. $\endgroup$
    – whuber
    Commented Jul 16, 2020 at 13:20
  • $\begingroup$ $det(R) >0$ is impossible for $|\rho_3|=1$ and $|\rho_2|=1$, why do you include them? Sylvester's criterion implies $|\rho_3|<1$, so it's impossible to generate a positive definite matrix $R$ with $|\rho_3|=1$. $\endgroup$
    – user289381
    Commented Jul 16, 2020 at 13:25
  • $\begingroup$ @ping I use the word "and," which means all the criteria must simultaneously apply. $\endgroup$
    – whuber
    Commented Jul 16, 2020 at 13:26
  • 2
    $\begingroup$ @AJKOER In addition to the working, efficient code, along with an example of its output that I provided, what more do you need for a demonstration? And, in this context, could you explain what you mean by "accuracy"? $\endgroup$
    – whuber
    Commented Jul 23, 2020 at 15:25

The G-Wishart distribution (Letac & Massam, 2007) is a distribution on positive definite matrices with fixed zeros corresponding to the missing edges of a graph $\mathcal G$ with nodes the indices $(i,j)$ of the associated variates. It has a density of the same form as the Wishart distribution: $$p(\Sigma|\delta,\Xi)\propto|\Sigma|^{(\delta-2)/2}\exp\left\{-\frac{1}{2}\text{tr}(\Sigma^\text{T}\Xi)\right\}$$and enjoys the most useful property that the conditional distributions of the submatrices of $\Sigma$ associated with the cliques of the graph all are standard Wishart, which allows for a Gibbs sampling approach to its simulation.

This distribution is implemented in R via the function rgwish. The graph $\mathcal G$ is described by an adjacency upper-triangular matrix adj that is made of 0's and 1's, with 0's indicating the fixed zeroes of the matrix.

In the current question, this R function can be called until all constraints are satisfied. The matrix $\Xi$ (denoted D in rgwish) can be chosen towards favouring the constraints to be met.

  • $\begingroup$ If you need a more flexible approach to test, say correlation robustness, to outliers, my spreadsheet approach may be of value as outlined in my edited thread. $\endgroup$
    – AJKOER
    Commented Jul 15, 2020 at 21:02
  • 1
    $\begingroup$ I challenge the potential feasibility and accuracy of the Wishart approach. Please demonstrate it, even in a simple example. $\endgroup$
    – AJKOER
    Commented Jul 23, 2020 at 14:38

Another possible solution is using eigendecomposition. You can use the following code to randomly generate an orthonormal matrix

In R:

Q_matrix = function(d=2){
  A = matrix(runif(d^2), d, d)
  Q = gramSchmidt(A)$Q # get orthonormal matrix by Gram Schmidt procedure (s.t. Q%*%t(Q) = I)

Then specify the diagonal matrix containing all the eigenvalues, L, then the covariance matrix will be

L = diag(1:2)
Q = Q_matrix(2)
S = Q%*%L%*%t(Q)

Here is how I would proceed to develop random normal deviates (Y, X, Z) which I would use to compute corresponding covariance matrices.

In the special case of a bivariate multinormal distribution, there is a formula that states for variable Y, which assumes a value of y (so Y=y), there is a conditional univariate normal distribution for X with mean as a function of y and the correlation between X and Y, and a formula for the variance. Here is a reference: https://www2.stat.duke.edu/courses/Spring12/sta104.1/Lectures/Lec22.pdf .

$E[X|Y = y] = µ_X + σ_X ρ\frac{(y − µ_Y)}{σ_Y} $

$Var[X|Y = y] = σ_X (1 − ρ^2) $

Assuming that the largest correlation exist between Y and X and also between Y and Z, then similarly:

$E[Z|Y = y] = µ_Z + σ_Z ρ\frac{(y − µ_Y)}{σ_Y} $

$Var[Z|Y = y] = σ_Z (1 − ρ^2) $

where the supplied correlations are applicable to Y and X, and also Y and Z, respectively.

So, using a package normal random deviate function, you simulate Y and use the value Y=y, to derive conditional univariate normal random deviate function values for X and Z, respectively.

[EDIT] To address a comment, my suggested approach generates data from which to construct a sample random covariance matrix. I would note, to quote Wikipedia:

Another issue is the robustness to outliers, to which sample covariance matrices are highly sensitive.[2][3][4]

However, if you have simulated point data, one may be able to address outliers (another topic) and perhaps relatedly mitigate this issue with your sample population of generated random covariance matrices.

[EDIT EDIT] Using a freely available google spreadsheet, I have a simulated 100 normal, and two sets (of 100) conditional normal random deviates (namely, X, Y and Z) with targeted parameters. Upon valuing the covariance matrix (and the associated correlation matrix), one can build a database of random simulated covariance matrices. Further, I can control the generating process to create a mixture of outliers, to investigate robustness estimates for covariance and correlations.

Here is a rough extract from the spreadsheet, please excuse the formating which is lost on pasting from the spreadsheet into a document text file.


Mean 0.25 0.1 0.2

SE 0.3 0.2 0.1

CORREL 1.0 0.7 0.5

OBS Mean 0.2177 0.0849 0.2028

OBS SE 0.3377 0.1576 0.0899

OBS VAR 0.1141 0.0081

OBS Correl 1 0.800 0.432

Example of underlying cell formulas employed:

${\text=NORMINV(RAND(),B$3,B$4) , =NORMINV(RAND(),0,C$4*SQRT(1-$C$5*$C$5)+C13 , =C$3+C$4*C$5*($B14-$B$3)/$B$4 , =E$3+E$4*E$5*($B14-$B$3)/$B$4 }$

0.1129 0.0538 0.0174
0.0538 0.0400 0.0077
0.0174 0.0077 0.0080

1.0000 0.8003 0.5777
0.8003 1.0000 0.4356
0.5777 0.4356 1.0000

Here is an output example of ranged valued Random Data Matrices produced:

/////////////////////////// DATA RUN DATABASE ////////////

0.0973 0.0529 0.0190
0.0529 0.0485 0.0102
0.0190 0.0102 0.0091

1.0000 0.7700 0.6373
0.7700 1.0000 0.5000
0.6373 0.500 1.0000

So far, the random matrice set I have created look pretty rationale without major swings from the targetted parameters. Has anyone seen the output from the Wishart approach, and would it be at all useful in testing robustness of say correlation matrices with a select noise component (as I can readily accomplish)?

  • $\begingroup$ To be clear, if one simulates a theoretical regression relationship between Y and X, one has, per that simulation run, a sample covariance matrix. Repeat the exercise k times and arrive at a set of k random covariance matrices, which should, at some point, meaningfully converge to a mean matrix and with an associated variability. The latter can be performed in a freely available google spreadsheet and can ID a matrix inverse approaching a singularity. In other words, my approach is not only intuitive but also transparent and robust approach to a sample Wishart distribution . $\endgroup$
    – AJKOER
    Commented Apr 12, 2020 at 10:48
  • $\begingroup$ Note, the problem cites correlations approaching 0.9, do you really expect no numerical analysis outliers to be created to distort the analysis? Construct a theoretical model and do test runs in whatever path you deem acceptable. $\endgroup$
    – AJKOER
    Commented Apr 12, 2020 at 11:06
  • $\begingroup$ I thanks those who challenged the feasibility of my flexible hands-on approach to generating random matrices capable of addressing my questions on, for example, the robustness of correlation matrices and/or covariance matrix by examining simulation runs with select outlier introduction (by, for example, creating a mixture of distributions in the generating step to introduce a noise model). $\endgroup$
    – AJKOER
    Commented Jul 15, 2020 at 20:48
  • $\begingroup$ So, I have demonstrated theoretical transparency, accessibility of implementation, and feasible results. Yet, the theoretical claim of a Wishart approach is unchallenged, excuse me, I hereby challenge its feasibility and accuracy. $\endgroup$
    – AJKOER
    Commented Jul 23, 2020 at 14:34
  • $\begingroup$ It's not clear what you mean by "feasible results." Your examples don't meet the constraints of the question on the correlation coefficients -- what else is there to say? $\endgroup$
    – whuber
    Commented Jul 23, 2020 at 15:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.