How to generate a large full-rank random correlation matrix with some strong correlations present? I would like to generate a random correlation matrix $\mathbf C$ of $n \times n$ size such that there are some moderately strong correlations present:


*

*square real symmetric matrix of $n \times n$ size, with e.g. $n=100$;

*positive-definite, i.e. with all eigenvalues real and positive;

*full rank;

*all diagonal elements equal to $1$;

*off-diagonal elements should be reasonably uniformly distributed on $(-1, 1)$. Exact distribution does not matter, but I would like to have some moderately large amount (e.g. $10\%$) of moderately large values (e.g. with absolute value of $0.5$ or higher). Basically I want to make sure that $\mathbf C$ is not almost diagonal with all off-diagonal elements $\approx 0$.


Is there a simple way to do it?
The purpose is to use such random matrices to benchmark some algorithms working with correlation (or covariance) matrices.

Methods that do not work
Here are some ways to generate random correlation matrices that I know of, but that do not work for me here:


*

*Generate random $\mathbf X$ of $s \times n$ size, center, standardize and form the correlation matrix $\mathbf C=\frac{1}{s-1}\mathbf X^\top \mathbf X$. If $s>n$, this will generally result in all off-diagonal correlations being around $0$. If $s\ll n$, some correlations will be strong, but $\mathbf C$ will not be full rank.

*Generate random positive definite matrix $\mathbf B$ in one of the following ways:


*

*Generate random square $\mathbf A$ and make symmetric positive definite $\mathbf B = \mathbf A \mathbf A^\top$.

*Generate random square $\mathbf A$, make symmetric $\mathbf E = \mathbf A + \mathbf A^\top$, and make it positive definite by performing eigen-decomposition $\mathbf E = \mathbf U \mathbf S \mathbf U^\top$ and setting all negative eigenvalues to zero: $\mathbf B = \mathbf U \:\mathrm{max}\{\mathbf S, \mathbf 0\} \:\mathbf U^\top$. NB: this will result in a rank-deficient matrix.

*Generate random orthogonal $\mathbf Q$ (e.g. by generating random square $\mathbf A$ and doing its QR decomposition, or via Gram-Schmidt process) and random diagonal $\mathbf D$ with all positive elements; form $\mathbf B = \mathbf Q \mathbf D \mathbf Q^\top$.
Obtained matrix $\mathbf B$ can be easily normalized to have all ones on the diagonal: $\mathbf C = \mathbf D^{-1/2}\mathbf B \mathbf D^{-1/2}$, where $\mathbf D = \mathrm{diag}\:\mathbf B$ is the diagonal matrix with the same diagonal as $\mathbf B$. All three ways listed above to generate $\mathbf B$ result in $\mathbf C$ having off-diagonal elements close $0$.

Update: Older threads
After posting my question, I found two almost duplicates in the past:


*

*How to generate random correlation matrix that has approximately normally distributed off-diagonal entries with given standard deviation?

*How to efficiently generate random positive-semidefinite correlation matrices?
Unfortunately, none of these threads contained a satisfactory answer (until now :)
 A: A simple thing but maybe will work for benchmark purposes: took your 2. and injected some correlations into starting matrix. Distribution is somewhat uniform, and changing $a$ you can get concentration near 1 and -1 or near 0. 
import numpy as np
from random import choice
import matplotlib.pyplot as plt

n = 100
a = 2

A = np.matrix([np.random.randn(n) + np.random.randn(1)*a for i in range(n)])
A = A*np.transpose(A)
D_half = np.diag(np.diag(A)**(-0.5))
C = D_half*A*D_half

vals = list(np.array(C.ravel())[0])
plt.hist(vals, range=(-1,1))
plt.show()
plt.imshow(C, interpolation=None)
plt.show()



A: Interesting question (as always!).  How about finding a set of example matrices that exhibit the properties you desire, and then take convex combinations thereof, since if $A$ and $B$ are positive definite, then so is $\lambda A + (1-\lambda)B$.  As a bonus, no rescaling of the diagonals will be necessary, by the convexity of the operation.  By adjusting the $\lambda$ to being more concentrated towards 0 and 1 versus uniformly distributed, you could concentrate the samples on the edges of the polytope, or the interior.  (You could use a beta/Dirichlet distribution to control the concentration vs uniformity).
For example, you could let $A$ to be component-symmetric, and $B$ be toeplitz.  Of course, you can always add another class $C$, and take $\lambda_A A + \lambda_B B + \lambda_C C$ such that $\sum \lambda = 1$ and $\lambda \geq 0$, and so on.
A: Other answers came up with nice tricks to solve my problem in various ways. However, I found a principled approach that I think has a large advantage of being conceptually very clear and easy to adjust.
In this thread: How to efficiently generate random positive-semidefinite correlation matrices? -- I described and provided the code for two efficient algorithms of generating random correlation matrices. Both come from a paper by Lewandowski, Kurowicka, and Joe (2009), that @ssdecontrol referred to in the comments above (thanks a lot!).
Please see my answer there for a lot of figures, explanations, and matlab code. The so called "vine" method allows to generate random correlation matrices with any distribution of partial correlations and can be used to generate correlation matrices with large off-diagonal values. Here is the example figure from that thread:

The only thing that changes between subplots, is one parameter that controls how much the distribution of partial correlations is concentrated around $\pm 1$.
I copy my code to generate these matrices here as well, to show that it is not longer than the other methods suggested here. Please see my linked answer for some explanations. The values of betaparam for the figure above were ${50,20,10,5,2,1}$ (and dimensionality d was $100$). 
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

Update: eigenvalues
@psarka asks about the eigenvalues of these matrices. On the figure below I plot the eigenvalue spectra of the same six correlation matrices as above. Notice that they decrease gradually; in contrast, the method suggested by @psarka generally results in a correlation matrix with one large eigenvalue, but the rest being pretty uniform.


Update. Really simple method: several factors
Similar to what @ttnphns wrote in the comments above and @GottfriedHelms in his answer, one very simple way to achieve my goal is to randomly generate several ($k<n$) factor loadings $\mathbf W$ (random matrix of $k \times n$ 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 (as described in my question). This is very simple and does the trick. Here are some example correlation matrices for $k={100, 50, 20, 10, 5, 1}$:

The only downside is that the resulting matrix will have $k$ large eigenvalues and then a sudden drop, as opposed to a nice decay shown above with the vine method. Here are the corresponding spectra:

Here is the code:
d = 100;    %// number of dimensions
k = 5;      %// number of factors

W = randn(d,k);
S = W*W' + diag(rand(1,d));
S = diag(1./sqrt(diag(S))) * S * diag(1./sqrt(diag(S)));

