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

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.

• Use Cholesky factorisation Apr 12 at 5:23

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.

• How can I construct back my matrix $M$? It seems that for $n=3$ in your code my matrix $M$ does not exist. Sorry if I misunderstood something.
– Matt
Apr 12 at 6:07
• If at the end you want $4\times4$-matrices $M$ correlated by a $4\times4\times4\times4$-tensor of elements $\tau_{ijkl}$, then set $n=16$, flatten the correlation tensor out into a $16\times16$-matrix T, and generate random 16-vectors. Then, take each 16-vector and reshape it into a $4\times4$ matrix. Apr 12 at 6:52