I want to generate $c$ correlated standard normals where the correlation matrix M is something like \begin{matrix} 1 & \rho & \rho \\ \rho & 1 & \rho \\ \rho & \rho & 1 \\ .... \end{matrix}
so ones along the diagonal but otherwise all variables have the same correlation $\rho$ where $\rho \geq 0$. What's the pattern in generating this using linear combinations of the independent standard normal variables?
Suppose that your covariance matrix $M$ has size $d\times d$. Then the gaussian distribution $\mathcal{N}(0,M)$ can be realized by the $d$-dimensional random vector $SZ$, where
In general, there are many ways of choosing such $S$ and often Cholesky decomposition is adopted. In our case, however, we can simply pick $S = M^{1/2}$ where $M^{1/2}$ is the unique positive-definite square root of $M$ as it turns out to be easily computed.
To compute $M^{1/2}$, we pick the column vector $u = \frac{1}{\sqrt{d}}(1, \cdots, 1)^{\mathsf{T}}$. Then $duu^{\mathsf{T}}$ is the $d\times d$ matrix consisting of only $1$'s and hence $M$ can be written as
$$M = (1-\rho)I_d + duu^{\mathsf{T}} = (1 - \rho)(I_d - uu^{\mathsf{T}}) + (1 + (d-1)\rho) uu^{\mathsf{T}}. \tag{*}$$
Notice that both $I_d - uu^{\mathsf{T}}$ and $uu^{\mathsf{T}}$ are orthogonal projections such that
$$(I_d - uu^{\mathsf{T}})uu^{\mathsf{T}} = uu^{\mathsf{T}}(I_d - uu^{\mathsf{T}}) = 0.$$
That is, $\text{(*)}$ is already a spectral decomposition of $M$. From this, we easily compute $M^{1/2}$ as
$$ M^{1/2} = \sqrt{\smash[b]{1-\rho}} \, (I_d - uu^{\mathsf{T}}) + \sqrt{\smash[b]{1 + (d-1)\rho}} \, uu^{\mathsf{T}}. $$
So if we put $\beta = \frac{1}{d}\left( \sqrt{\smash[b]{1 + (d-1)\rho}} - \sqrt{\smash[b]{1-\rho}} \right)$ and $\alpha = \sqrt{\smash[b]{1-\rho}} + \beta$, then
$$ M^{1/2} = \begin{pmatrix} \alpha & \beta & \beta & \cdots & \beta \\ \beta & \alpha & \beta & \cdots & \beta \\ \beta & \beta & \alpha & \cdots & \beta \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \beta & \beta & \beta & \cdots & \alpha \end{pmatrix}. $$
