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.


*

*Could you suggest any efficient solutions? If you are aware of any examples in Matlab, I would be very thankful.

*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?

 A: 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.
A: 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.
A: 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.
A: 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:


*

*Rotations.

*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: 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).
A: 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)

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

A: 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.
