How to efficiently generate random positive-semidefinite correlation matrices?

I would like to be able to efficiently generate positive-semidefinite (PSD) correlation matrices. My method slows down dramatically as I increase the size of matrices to be generated.

1. Could you suggest any efficient solutions? If you are aware of any examples in Matlab, I would be very thankful.
2. When generating a PSD correlation matrix how would you pick the parameters to describe matrices to be generated? An average correlation, standard deviation of correlations, eigenvalues?

You can do it backward: every matrix $C \in \mathbb{R}_{++}^p$ (the set of all symmetric $p \times p$ PSD matrices) can be decomposed as

$C=O^{T}DO$ where $O$ is an orthonormal matrix

To get $O$, first generate a random basis $(v_1,...,v_p)$ (where $v_i$ are random vectors, typically in $(-1,1)$). From there, use the Gram-Schmidt orthogonalization process to get $(u_1,....,u_p)=O$

$R$ has a number of packages that can do the G-S orthogonalization of a random basis efficiently, that is even for large dimensions, for example the 'far' package. Although you will find the G-S algorithm on wiki, it's probably better not to re-invent the wheel and go for a matlab implementation (one surely exists, i just can't recommend any).

Finally, $D$ is a diagonal matrices whose elements are all positive (this is, again, easy to generate: generate $p$ random numbers, square them, sort them and place them unto the diagonal of a identity $p$ by $p$ matrix).

• (1) Note that the resulting $C$ will not be a correlation matrix (as requested by the OP), because it will not have ones on the diagonal. Of course it can be rescaled to have ones on the diagonal, by setting it to $E^{-1/2}CE^{-1/2}$, where $E$ is a diagonal matrix with the same diagonal as $C$. (2) If I am not mistaken, this will result in correlation matrices where all the off-diagonal elements are concentrated around $0$, so there is no flexibility that the OP was looking for (OP wanted to be able to set "an average correlation, standard deviation of correlations, eigenvalues") – amoeba Nov 19 '14 at 10:48
• @amoeba: I will address (2) since, as you point out, the solution to (1) is trivial. The one number characterization of the 'shape' (the relationship between the in and out of diagonal elements) of a PSD matrix (and hence a covariance and also a correlation matrix) is its condition number. And, the method described above allows full control over it. The 'concentration of the off diagonal elements around 0' is not a feature of the method used to generate PSD matrices but, rather, a consequence of the constraints necessary to ensure that the matrix is PSD and the fact that $p$ is large. – user603 Nov 19 '14 at 12:57
• Are you saying that all large PSD matrices have off-diagonal elements close to zero? I disagree, it is not so. Check my answer here for some examples: How to generate random correlation matrix that has approximately normally distributed off-diagonal entries with given standard deviation? But one can directly see that it is not the case, because a square matrix having all ones on the diagonal and a fixed value $\rho$ everywhere off-diagonal is PSD and $\rho$ can be arbitrarily large (but of course below $1$). – amoeba Nov 19 '14 at 13:14
• @amoeba: then I was wrong in assuming that by necessity the off diagonal of large correlation matrices when they are allowed to be both positive as well as negative are close to 0. Thanks for the enlightening example. – user603 Nov 19 '14 at 13:56
• I read a very nice paper on generating random correlation matrices and provided my own answer here (as well as another answer in that linked thread). I think you might find it interesting. – amoeba Nov 21 '14 at 23:18

A paper Generating random correlation matrices based on vines and extended onion method by Lewandowski, Kurowicka, and Joe (LKJ), 2009, provides a unified treatment and exposition of the two efficient methods of generating random correlation matrices. Both methods allow to generate matrices from a uniform distribution in a certain precise sense defined below, are simple to implement, fast, and have an added advantage of having amusing names.

A real symmetric matrix of $d \times d$ size with ones on the diagonal has $d(d-1)/2$ unique off-diagonal elements and so can be parametrized as a point in $\mathbb R^{d(d-1)/2}$. Each point in this space corresponds to a symmetric matrix, but not all of them are positive-definite (as correlation matrices have to be). Correlation matrices therefore form a subset of $\mathbb R^{d(d-1)/2}$ (actually a connected convex subset), and both methods can generate points from a uniform distribution over this subset.

I will provide my own MATLAB implementation of each method and illustrate them with $d=100$.

Onion method

The onion method comes from another paper (ref #3 in LKJ) and owns its name to the fact the correlation matrices are generated starting with $1\times 1$ matrix and growing it column by column and row by row. Resulting distribution is uniform. I don't really understand the math behind the method (and prefer the second method anyway), but here is the result:

Here and below the title of each subplot shows the smallest and the largest eigenvalues, and the determinant (product of all eigenvalues). Here is the code:

%// ONION METHOD to generate random correlation matrices distributed randomly
function S = onion(d)
S = 1;
for k = 2:d
y = betarnd((k-1)/2, (d-k)/2); %// sampling from beta distribution
r = sqrt(y);
theta = randn(k-1,1);
theta = theta/norm(theta);
w = r*theta;
[U,E] = eig(S);
R = U*E.^(1/2)*U';             %// R is a square root of S
q = R*w;
S = [S q; q' 1];               %// increasing the matrix size
end
end


Extended onion method

LKJ modify this method slightly, in order to be able to sample correlation matrices $\mathbf C$ from a distribution proportional to $[\mathrm{det}\:\mathbf C]^{\eta-1}$. The larger the $\eta$, the larger will be the determinant, meaning that generated correlation matrices will more and more approach the identity matrix. The value $\eta=1$ corresponds to uniform distribution. On the figure below the matrices are generated with $\eta={1, 10, 100, 1000, 10\:000, 100\:000}$.

For some reason to get the determinant of the same order of magnitude as in the vanilla onion method, I need to put $\eta=0$ and not $\eta=1$ (as claimed by LKJ). Not sure where the mistake is.

%// EXTENDED ONION METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = extendedOnion(d, eta)
beta = eta + (d-2)/2;
u = betarnd(beta, beta);
r12 = 2*u - 1;
S = [1 r12; r12 1];

for k = 3:d
beta = beta - 1/2;
y = betarnd((k-1)/2, beta);
r = sqrt(y);
theta = randn(k-1,1);
theta = theta/norm(theta);
w = r*theta;
[U,E] = eig(S);
R = U*E.^(1/2)*U';
q = R*w;
S = [S q; q' 1];
end
end


Vine method

Vine method was originally suggested by Joe (J in LKJ) and improved by LKJ. I like it more, because it is conceptually easier and also easier to modify. The idea is to generate $d(d-1)/2$ partial correlations (they are independent and can have any values from $[-1, 1]$ without any constraints) and then convert them into raw correlations via a recursive formula. It is convenient to organize the computation in a certain order, and this graph is known as "vine". Importantly, if partial correlations are sampled from particular beta distributions (different for different cells in the matrix), then the resulting matrix will be distributed uniformly. Here again, LKJ introduce an additional parameter $\eta$ to sample from a distribution proportional to $[\mathrm{det}\:\mathbf C]^{\eta-1}$. The result is identical to the extended onion:

%// VINE METHOD to generate random correlation matrices
%// distributed ~ det(S)^eta [or maybe det(S)^(eta-1), not sure]
function S = vine(d, eta)
beta = eta + (d-1)/2;
P = zeros(d);           %// storing partial correlations
S = eye(d);

for k = 1:d-1
beta = beta - 1/2;
for i = k+1:d
P(k,i) = betarnd(beta,beta); %// sampling from beta
P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
p = P(k,i);
for l = (k-1):-1:1 %// converting partial correlation to raw correlation
p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
end
S(k,i) = p;
S(i,k) = p;
end
end
end


Vine method with manual sampling of partial correlations

As one can see above, uniform distribution results in almost-diagonal correlation matrices. But one can easily modify the vine method to have stronger correlations (this is not described in the LKJ paper, but is straightforward): for this one should sample partial correlations from a distribution concentrated around $\pm 1$. Below I sample them from beta distribution (rescaled from $[0,1]$ to $[-1, 1]$) with $\alpha=\beta={50, 20, 10, 5, 2, 1}$. The smaller the parameters of the beta distribution, the more it is concentrated near the edges.

Note that in this case the distribution is not guaranteed to be permutation invariant, so I additionally randomly permute rows and columns after generation.

%// VINE METHOD to generate random correlation matrices
%// with all partial correlations distributed ~ beta(betaparam,betaparam)
%// rescaled to [-1, 1]
function S = vineBeta(d, betaparam)
P = zeros(d);           %// storing partial correlations
S = eye(d);

for k = 1:d-1
for i = k+1:d
P(k,i) = betarnd(betaparam,betaparam); %// sampling from beta
P(k,i) = (P(k,i)-0.5)*2;     %// linearly shifting to [-1, 1]
p = P(k,i);
for l = (k-1):-1:1 %// converting partial correlation to raw correlation
p = p * sqrt((1-P(l,i)^2)*(1-P(l,k)^2)) + P(l,i)*P(l,k);
end
S(k,i) = p;
S(i,k) = p;
end
end

%// permuting the variables to make the distribution permutation-invariant
permutation = randperm(d);
S = S(permutation, permutation);
end


Here is how the histograms of the off-diagonal elements look for the matrices above (variance of the distribution monotonically increases):

Update: using random factors

One really simple method of generating random correlation matrices with some strong correlations was used in the answer by @shabbychef, and I would like to illustrate it here as well. The idea is to randomly generate several ($k<d$) factor loadings $\mathbf W$ (random matrix of $k \times d$ size), form the covariance matrix $\mathbf W \mathbf W^\top$ (which of course will not be full rank) and add to it a random diagonal matrix $\mathbf D$ with positive elements to make $\mathbf B = \mathbf W \mathbf W^\top + \mathbf D$ full rank. The resulting covariance matrix can be normalized to become a correlation matrix, by letting $\mathbf C = \mathbf E^{-1/2}\mathbf B \mathbf E^{-1/2}$, where $\mathbf E$ is a diagonal matrix with the same diagonal as $\mathbf B$. This is very simple and does the trick. Here are some example correlation matrices for $k={100, 50, 20, 10, 5, 1}$:

And the code:

%// FACTOR method
function S = factor(d,k)
W = randn(d,k);
S = W*W' + diag(rand(1,d));
S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));
end


Here is the wrapping code used to generate the figures:

d = 100; %// size of the correlation matrix

figure('Position', [100 100 1100 600])
for repetition = 1:6
S = onion(d);

%// etas = [1 10 100 1000 1e+4 1e+5];
%// S = extendedOnion(d, etas(repetition));

%// S = vine(d, etas(repetition));

%// betaparams = [50 20 10 5 2 1];
%// S = vineBeta(d, betaparams(repetition));

subplot(2,3,repetition)

%// use this to plot colormaps of S
imagesc(S, [-1 1])
axis square
title(['Eigs: ' num2str(min(eig(S)),2) '...' num2str(max(eig(S)),2) ', det=' num2str(det(S),2)])

%// use this to plot histograms of the off-diagonal elements
%// offd = S(logical(ones(size(S))-eye(size(S))));
%// hist(offd)
%// xlim([-1 1])
end

• This is a fantastic rundown, I'm glad I said something! – shadowtalker Nov 22 '14 at 11:05
• When I translated the matlab code for the vine based correlation matrix to R and tested it, the density of correlations in column 1 was always different to later columns. It may be that I translated something incorrectly, but perhaps this note helps someone. – Charlie Nov 1 '16 at 22:15
• For R users, the function rcorrmatrix in package clusterGeneration (written by W Qui and H. Joe) implements the vine method. – Robert Montgomery Jul 12 '18 at 16:24

An even simpler characterization is that for real matrix $A$, $A^TA$ is positive semidefinite. To see why this is the case, one only has to prove that $y^T (A^TA) y \ge 0$ for all vectors $y$ (of the right size, of course). This is trivial: $y^T (A^TA) y = (Ay)^T Ay = ||Ay||$ which is nonnegative. So in Matlab, simply try

A = randn(m,n);   %here n is the desired size of the final matrix, and m > n
X = A' * A;


Depending on the application, this may not give you the distribution of eigenvalues you want; Kwak's answer is much better in that regard. The eigenvalues of X produced by this code snippet should follow the Marchenko-Pastur distribution.

For simulating the correlation matrices of stocks, say, you may want a slightly different approach:

k = 7;      % # of latent dimensions;
n = 100;    % # of stocks;
A = 0.01 * randn(k,n);  % 'hedgeable risk'
D = diag(0.001 * randn(n,1));   % 'idiosyncratic risk'
X = A'*A + D;
ascii_hist(eig(X));    % this is my own function, you do a hist(eig(X));
-Inf <= x <  -0.001 : **************** (17)
-0.001 <= x <   0.001 : ************************************************** (53)
0.001 <= x <   0.002 : ******************** (21)
0.002 <= x <   0.004 : ** (2)
0.004 <= x <   0.005 :  (0)
0.005 <= x <   0.007 : * (1)
0.007 <= x <   0.008 : * (1)
0.008 <= x <   0.009 : *** (3)
0.009 <= x <   0.011 : * (1)
0.011 <= x <     Inf : * (1)

• would you be willing to share your ascii_hist function by any chance? – btown Jan 1 '11 at 0:02
• @btown the margin is too small to contain it! – shabbychef Jan 2 '11 at 5:34
• There is a typo in $y^T (A^TA) y = (Ay)^T Ay = ||Ay||$ - it's missing its final square! – Silverfish Dec 26 '14 at 21:06

As a variation on kwak's answer: generate a diagonal matrix $\mathbf{D}$ with random nonnegative eigenvalues from a distribution of your choice, and then perform a similarity transformation $\mathbf{A}=\mathbf{Q}\mathbf{D}\mathbf{Q}^T$ with $\mathbf{Q}$ a Haar-distributed pseudorandom orthogonal matrix.

• M. : Nice reference: This appears to be the most efficient solution (asymptotically). – whuber Sep 17 '10 at 14:46
• @whuber: Heh, I picked it up from Golub and Van Loan (of course); I use this all the time to help in generating test matrices for stress-testing eigenvalue/singular value routines. As can be seen from the paper, it's essentially equivalent to QR-decomposing a random matrix like what kwak suggested, except that it is done more efficiently. There's a MATLAB implementation of this in Higham's Text Matrix Toolbox, BTW. – J. M. is not a statistician Sep 17 '10 at 19:14
• M.:> Thanks for the matlab implementation. Would you by any chance know of a Haar pseudo-random matrix generator in R? – user603 Sep 17 '10 at 19:58
• @kwak: No idea, but if there's no implementation yet, it shouldn't be too hard to translate the MATLAB code to R (I can try to whip one up if there really isn't any); the only prerequisite is a decent generator for pseudorandom normal variates, which I'm sure R has. – J. M. is not a statistician Sep 17 '10 at 20:06
• M.:> yes i'll probably translate it my self. Thanks for the links, Best. – user603 Sep 17 '10 at 20:13

You haven't specified a distribution for the matrices. Two common ones are the Wishart and inverse Wishart distributions. The Bartlett decomposition gives a Cholesky factorisation of a random Wishart matrix (which can also be efficiently solved to obtain a random inverse Wishart matrix).

In fact, the Cholesky space is a convenient way to generate other types of random PSD matrices, as you only have to ensure that the diagonal is non-negative.

• > Not random: two matrices generated from the same Whishard will not be independant from one another. If you plan to change the Whishart at each generation, then how do you intend to generate those Whishart in the first place? – user603 Sep 24 '10 at 9:32
• @kwak: I don't understand your question: the Bartlett decomposition will give independent draws from the same Wishart distribution. – Simon Byrne Sep 24 '10 at 11:03
• > Let me rephrase this, where do you get the scale matrix of your whishart distribution from ? – user603 Sep 24 '10 at 12:16
• @kwak: it's a parameter of the distribution, and so is fixed. You select it at the start, based on the desired characteristics of your distribution (such as the mean). – Simon Byrne Sep 24 '10 at 12:40

The simplest method is the one above, which is a simulation of a random dataset and the computation of the Gramian. A word of caution: The resulting matrix will not be uniformly random, in that its decomposition, say $U^TSU$ will have rotations not distributed according to the Haar Measure. If you want to have "uniformly distributed" PSD matrices then you can use any of the approaches described here.

• If the entries are generated from a Normal distribution rather than a uniform, the decomposition you mention ought to be SO(n) invariant (and therefore equidistributed relative to the Haar measure). – whuber Sep 17 '10 at 11:47
• Interesting. Can you provide a reference for this? – gappy Sep 17 '10 at 14:39
• > the problem with this method is that you can't control the ratio of smallest to largest eigenvalue (and i think that as the size of your randomly generated dataset goes to infinity, this ratio will converge to 1). – user603 Sep 24 '10 at 9:28

If you'd like to have more control over your generated symmetric PSD matrix, e.g. generate a synthetic validation dataset, you have a number of parameters available. A symmetric PSD matrix corresponds to a hyper-ellipse in the N-dimensional space, with all the related degrees of freedom:

1. Rotations.
2. Lengths of axes.

So, for a 2-dimensional matrix (i.e. 2d ellipse), you'll have 1 rotation + 2 axes = 3 parameters.

If rotations bring to mind Orthogonal matrices, its a correct train of though, since the construction is again $\Sigma=ODO^T$, with $\Sigma$ being the produced Sym.PSD matrix, $O$ the rotation matrix (which is orthogonal), and $D$ the diagonal matrix, whose diagonal elements will control the length of the axes of the ellipse.

The following Matlab code plots 16 2dimensional Gaussian-distributed datasets based on $\Sigma$, with an increasing angle. The code for random generation of parameters is in comments.

figure;
mu = [0,0];
for i=1:16
subplot(4,4,i)
theta = (i/16)*2*pi;   % theta = rand*2*pi;
U=[cos(theta), -sin(theta); sin(theta) cos(theta)];
% The diagonal's elements control the lengths of the axes
D = [10, 0; 0, 1]; % D = diag(rand(2,1));
sigma = U*D*U';
data = mvnrnd(mu,sigma,1000);
plot(data(:,1),data(:,2),'+'); axis([-6 6 -6 6]); hold on;
end


For more dimensions, the Diagonal matrix is straight-forward (as above), and the $U$ should derive from multiplication of the rotation matrices.

A cheap and cheerful approach I've used for testing is to generate m N(0,1) n-vectors V[k] and then use P = d*I + Sum{ V[k]*V[k]'} as an nxn psd matrix. With m < n this will be singular for d=0, and for small d will have high condition number.

• > the problem with this method is that you can't control the ratio of smallest to largest eigenvalue (and i think that as the size of your randomly generated dataset goes to infinity, this ratio will converge to 1). – user603 Sep 24 '10 at 9:28
• > besides, the method is not very efficient (from a computationnal point of view) – user603 Sep 24 '10 at 9:36
• Your "random matrix" is a specially structured one called a "diagonal plus rank-1 matrix" (DR1 matrix), so not really a good representative random matrix. – J. M. is not a statistician Sep 24 '10 at 10:08