# Create a $N\times N$ random matrix with correlated off-diagonal elements

I am working on spectral distribution of random matrices. I need to use simulations and I have been switching between MATLAB and Mathematica. I wish to perform all my computations on Mathematica.

In order to do so, I would like to create a random matrix $$\mathbf{A}$$ of size $$N\times N$$ such that all $$A_{ij}$$ are identically distributed with mean $$\mathbb{E}[A_{ij}]=0$$ and variance $$\frac{1}{N}$$, and $$A_{ij}$$ is correlated with $$A_{ji}$$ i.e $$\mathbb{E}[A_{ij}A_{ji}]\cdot N=\tau$$. (I have to scale the variance by $$N^{-1}$$ otherwise the eigenvalues will be spread on the entire complex plane).

This can be done in a few lines with MATLAB but I am encountering a lot of difficulties implementing it on Mathematica.

Any help or advice is highly appreciated, thank you.

To generate the correlated pairs you could use the BinormalDistribution function:

dist = BinormalDistribution[{0, 0}, {1/Sqrt[n], 1/Sqrt[n]}, τ];
Expectation[n x1 x2 , {x1, x2} \[Distributed] dist]
(* τ *)


To put this all together to obtain a realization of the random matrix $$\boldsymbol{A}$$, the following brute force approach could work. (Only $$n(n+1)/2$$ pairs are necessary to generate.)

matrix[n_, τ_] := Module[{x},
(* Random sample of correlated pairs *)
x = RandomVariate[BinormalDistribution[{0, 0}, {1/Sqrt[n], 1/Sqrt[n]}, τ], n (n + 1)/2];

(* Construct n x n matrix from the random sample *)
k = 0;
a = ConstantArray[0, {n, n}];
Do[
Do[
k = k + 1;
If[i == j, a[[i, j]] = x[[k, 1]],
a[[i, j]] = x[[k, 1]]; a[[j, i]] = x[[k, 2]]],
{j, i, n}],
{i, 1, n}];

(* Return matrix *)
a
]

SeedRandom[12345];
matrix[5, 1/2] // MatrixForm


$$\left( \begin{array}{ccccc} -0.442184 & 0.805261 & -0.283361 & 0.132974 & -0.753305 \\ -0.153694 & 0.178089 & -0.0917344 & -0.336303 & 0.241191 \\ -0.49796 & 0.119943 & 0.188825 & -0.363142 & 0.504677 \\ -0.443263 & 0.0284703 & -0.316425 & -0.05942 & -0.488706 \\ -0.940153 & 0.586686 & -0.211216 & -0.272132 & -0.351116 \\ \end{array} \right)$$

To generate two Gaussian random numbers that both have zero mean and unit variance, and that have a covariance of $$\tau\in[-1,1]$$, you can do

twoCorrelatedGaussianRandoms[τ_] :=
({{Sqrt[1+τ]+Sqrt[1-τ], Sqrt[1+τ]-Sqrt[1-τ]},
{Sqrt[1+τ]-Sqrt[1-τ], Sqrt[1+τ]+Sqrt[1-τ]}}/2) . RandomVariate[NormalDistribution[], 2]


Try it out by generating a million pairs with $$\tau=0.2$$:

P = Table[twoCorrelatedGaussianRandoms[0.2], 10^6];
P // Mean
(*    {-0.00129118, -0.000991929}    *)
P // Covariance
(*    {{1.00194, 0.200685}, {0.200685, 1.00014}}    *)


You can construct your matrix from this.