How to generate a random matrix with arbitrary correlation between elements?

Question

I would like to find a smart way to generate a $N\times N$ random matrix $M$ with arbitrary correlation: \begin{equation} \boxed{\langle M_{ij}M_{kl}\rangle=\tau_{ijkl}} \end{equation} Where the mean and variance of the elements are given by: \begin{align} \langle M_{ij}\rangle&=0 \\ \langle M_{ij}^2\rangle&=\sigma^2 \end{align}

---

The case I am interested in is actually a sub-problem of this. I would like to generate a matrix whose elements follow a normal distribution of mean $0$ and variance $1/N$, and whose elements are correlated the following way: \begin{equation} \langle M_{ij}M_{ki}\rangle=\tau_{ijk} \end{equation} When $\tau_{ijk}=\delta_{jk}N^{-1}$ I recover a symmetric matrix.

Answer 1

Simplify $ij\to u$ and $kl\to v$: find $M_u$ such that $\langle M_u\rangle=0$ and $\langle M_u M_v\rangle=\tau_{u v}$. This notation makes the analysis a bit simpler.

We can achieve this effect by using the matrix-square-root of the matrix $T$ of elements $\tau_{u v}$. Example:

n = 3;
T = {{2, -1, 0.3}, {-1, 4, 1.3}, {0.3, 1.3, 2}};
SymmetricMatrixQ[T] && PositiveSemidefiniteMatrixQ[T]
(*    True    *)

Using the matrix square root MatrixPower[T, 1/2] to generate $10^5$ lists of random numbers and computing their covariance matrix:

V = RandomVariate[NormalDistribution[], {10^5, n}] . MatrixPower[T, 1/2];
Covariance[V]
(*    {{1.99821, -1.01062, 0.29987},
       {-1.01062, 4.03928, 1.28961},
       {0.29987, 1.28961, 1.99284}}    *)

We see that this covariance matrix matches the desired T.