Let $M$ be a $N\times N$ random symmetric matrix, which is drawn from a Gaussian distribution of mean $\mu=0$ and variance $\sigma^2=1/N$.

The first level correlation is defined by: $$\begin{equation} \rho(\lambda)=\left \langle \text{Tr}\delta(\lambda-M)\right \rangle\quad\quad\quad\quad \tag{1} \label{eqn1} \end{equation}$$ Which is nothing more than the distribution of eigenvalues: following my routine below I can compute this quantity. The brackets represents the average over the gaussian distribution.

RandomMatrix[n_] := 
  RandomMatrix[n] = 
   RandomVariate[NormalDistribution[0, 1/Sqrt[n]], {n, n}];
n := 1500
M := 1/Sqrt[2]*(RandomMatrix[n] + Transpose@RandomMatrix[n])
\[Rho] = Eigenvalues[M];

enter image description here

In the large limit $N\gg 1$ equation \eqref{eqn1} will describe a semi-circular law.


The two-level correlation function is defined by: $$\begin{equation} \rho^{(2)}(\lambda, \mu)=\left\langle\frac{1}{N} \operatorname{Tr} \delta(\lambda-M) \frac{1}{N} \operatorname{Tr} \delta(\mu-M)\right\rangle\quad\quad\quad\quad \tag{2} \label{eqn2} \end{equation}$$ I interpret this quantity is as how correlated two given eigenvalues will be. Following the previous routine, my $\rho$ will be a list of the $N$ eigenvalues of $M$. How can I numerically compute equation \eqref{eqn2}? Isn't it impossible to tell how correlated two elements are in one same list? If I generate several lists, it has to be for different matrices, which does not work anymore.

Given one matrix $M$, how can I compute the correlation between two eigenvalues?

Any help or remark is appreciated!

  • 1
    $\begingroup$ Your definition of the square random matrix gives the same random matrix each time—is that what you intend? Here you are looking at distributions for a single random matrix and not a sample of many random matrices of a given size. Also, you don’t need to use set-delayed in your definition for n. The set-delaleyed for M is the right idea of you want different random matrices each time, but won’t work with your definition of RandomMatrix. $\endgroup$ Jun 21, 2021 at 8:27
  • $\begingroup$ Yes it was intended. For two reasons: first because when I want to symmetrize M I need RandomMatrix[n] to stay the same when I call it twice: RandomMatrix[n]+Transpose@RandomMatrix[n]. The second reason was to emphasize the fact that when I speak of computing the correlation function, it is for a single matrix $M$ and not a collection of random matrices. For example when I computed the eigenvalues it was for one matrix, I didn't average it over many matrices. $\endgroup$
    – Matt
    Jun 21, 2021 at 9:45
  • $\begingroup$ However I understand from your comment that I should not use $n$ as the set-delayed, I will use a different parameter next time, thank you. $\endgroup$
    – Matt
    Jun 21, 2021 at 9:47
  • $\begingroup$ One point is hard to understand from your post, why don't you compute according to Eq.(1), but do instead a matrix diagonalization? $\endgroup$
    – yarchik
    Jun 21, 2021 at 10:08
  • $\begingroup$ One minor point - the diagonal elements of your RandomMatrix don't have the same distribution as the off-diagonal element (Variance[Diagonal[M]] / Variance[Flatten[M]] is about 2). The off-diagonals are the sum of 2 independent RVs, while the diagonal is the sum of the same RV. $\endgroup$
    – mikado
    Jun 21, 2021 at 10:58

1 Answer 1


The question is: How do you calculate the two-level correlation function for Gaussian Orthogonal Ensemble (GOE)? One way would be to create a 2d histogram. You create bins of some size $\Delta\lambda \times \Delta\mu$ and then for each bin add the number of eigenvalues in the interval $\lambda + \Delta\lambda$ to the number of eigenvalues in the interval $\mu + \Delta\mu$. The bins should be small enough to see the structure you are interested in, but large enough to capture enough eigenvalues to get decent numerics. I'm not sure it would be easy to get good numerics with just one matrix. For the self-averaging property of the GOE to give a smooth two-point correlation function, the matrix would have to be very big, I guess maybe $10^8$ by $10^8$, but it would take a long time and a lot of RAM to diagonalize such a large matrix. It would be better to do an ensemble average over smaller matrices.

I don't know what 2d histogram features are available in Mathematica. I looked into calculating a 2d histogram of the two point correlation function of the Gaussian Unitary Ensemble (GUE) using Matlab or Octave which have some built-in 2d histogram functions. But in the end I decided to use the formula for $\rho^{(2)}(\lambda,\mu)$ explicitly with gaussian-broadened delta functions. The widths of gaussian functions used to represent the delta functions correspond to the bin sizes in the histogram approach. Here is the resulting publication: https://link.springer.com/article/10.1007/s13538-020-00802-6

  • 2
    $\begingroup$ This does not provide an answer to the question. Once you have sufficient reputation you will be able to comment on any post; instead, provide answers that don't require clarification from the asker. - From Review $\endgroup$
    – corey979
    Apr 18, 2023 at 14:56
  • 1
    $\begingroup$ @corey979 The asker asked "How can I numerically compute equation (2)?" I answered that this could be done by creating a 2d histogram and I described a procedure for coding the 2d histogram in any language. Built-in Mathematica functions may produce an acceptable 2d histogram of Eq. (2) but it is more likely that the asker would have code the histogram himself. In my answer I also described an alternative manner of calculating Eq. (2) using gaussian broadened delta functions which also could be coded in any language. $\endgroup$ Apr 18, 2023 at 20:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.