How to efficiently generate random positive-semidefinite correlation matrices?

Question

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?

Answer 1

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).

Answer 2

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

Answer 3

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)

Answer 4

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.

Answer 5

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.

Answer 6

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.

Answer 7

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.

Answer 8

If you want to sample random correlation matrices from an empirical distribution, you can try to use generative models from machine learning to do that.

One such example is CorrGAN: Sampling Realistic Financial Correlation Matrices Using Generative Adversarial Networks

The basic idea is to fit generative adversarial networks (or variational autoencoders) to a set of empirical correlation matrices (having certain properties that are hard to capture/generate mathematically).

Then, you can use these models to sample as much random correlation matrices as you want that verify the properties of the empirical ones.

Depending on which networks you use, you may want to project the output (which may be not totally PSD) onto the set of correlation matrices (with the Higham projection for example).

Answer 9

you can create an arbitrary covariance matrix from the Wishart distribution using the stats function rWishart (included in base R) and then transform it to a correlation matrix

p <- 4
S <- drop(rWishart(1,p,diag(p)))
S * (diag(S)^(-1/2) %o% diag(S)^(-1/2))

The Wishart distribution takes two parameters: the degrees of freedom $\nu$ and the scale matrix $\Sigma$. To understand how these influence the result, condier $X_i \overset{iid}{\sim} N_p(0,\Sigma)$,then $$ S = \sum\limits_{i=1}^\nu X_iX_i' \sim W_p(\nu,\Sigma) $$ That is, S follows a Wishart distribution on $p\times p$ matrices with $\nu$ degrees of freedom and $\Sigma$ as scale matrix. As seen, if $\nu$ is larger than $p$, then $S$ will be positive definite a.s.(P). Furthermore, $$ E(S) = \nu \Sigma $$ which means that on average, the random draws S will approach a value proportional to the true population parameters $\Sigma$ in probability. This rate of convergence depends on the value of $\nu$, so if $p$ is very large, try to keep $\nu$ close to $p$, otherwise the random draws would be close to each other.

Answer 10

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.