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)));
The number of eigenvalues itself tell us nothing that we did not know (probably) already, namely the dimensions of the original matrix $A$ (or the dimensions of system of equations represented by $A$ depends how we see things).
The number of unique eigenvalues can informative (eg. the number of unique solutions to the system of equation defined by the matrix $A$), the number of non-zero eigenvalues can be informative (eg. the rank of the matrix $A$), the number of positive or negative eigenvalues can be informative (eg. is $A$ non-negative, negative, etc.), their type being real or complex can be informative (eg. is $A$ symmetric - the eigenvalues of a symmetric matrix are always real), their magnitude compared to each other can be informative ( eg. the condition number of $A$, compactness of $A$'s spectrum, etc.) but actually their number.. not much really.
The only things I would draw extra attention on is that the number of eigenvalues is NOT the cardinality of the spectrum, which is the number of unique eigenvalues. In general, recall that the fundamental theorem of algebra guarantees that every polynomial of degree one or more has a possibly complex root and we are just solving the characteristic polynomial of $A$ such that: $p_A(\lambda)=|A- \lambda I|=0$ to get the eigenvalues. The fact got $n$ of them (with multiplicity) doesn't say much.
As whuber commented, please note that a single value might occur as a root of the characteristic equation $\nu$ times in which case we say that value has (algebraic) multiplicity $\nu$. To that extent, when we deal with a $n \times n$ matrix we have $n$-th degree polynomial and as such the sum of the multiplicities of distinct eigenvalues is $n$.
This thread on Eigenvalues are unique? in Math.SE is also informative on the matter.
Best Answer
Just some hint to the first question, the remaining three are tractable, but require more context and answers to them are probably more verbose. As I commented, it is better to ask them in separate posts and show some of your own work.
You can write $\Delta$ more explicitly as $\Delta = \operatorname{diag}(a_1, \ldots, a_M)$ with $a_i > 0, i = 1, \ldots, M$. Denote the trace of $\Delta$ by $T = a_1 + \cdots + a_M$. It can be shown that $L_\Delta = \Delta^{1/2}C\Delta^{1/2}$, where \begin{align} C = I_{(M)} - \frac{1}{T}vv', \quad v = \begin{bmatrix} \sqrt{a_1} \\ \sqrt{a_2} \\ \vdots \\ \sqrt{a_M} \end{bmatrix}. \end{align}
Therefore to show $L_\Delta$ is PSD, it suffices to show $C$ is PSD. But $C$ is symmetric and idempotent, hence PSD. This completes the proof.