Eigenvalue Density via the Stieltjes Transform

These are my notes on how to compute the eigenvalue density of a random matrix using the Stieltes transform. They’re based on Chapter 2 of “A First Course in Random Matrix Theory,” and conversations with ChatGPT. I am new to this topic so these notes may contain errors.

We’re interested in the distribution of eigenvalues of a symmetric $N \times N$ random matrix $X_N$. We’ll ultimately be interested in the limit of infinitely large $N$, but we’ll start with the finite case.

There are lots of ways we can encode the eigenvalues. One natural way is to encode them as a probability distribution. Since we (for now) only have a discrete set of eigenvalues, we can put a delta function at each eigenvalue, so $$ \rho_N(\lambda) = {1 \over N} \sum_{i=1}^N \delta(\lambda – \lambda_i).$$

We will want to bring in tools from complex analysis, so another way to encode the eigenvalues is to create an “indicator” function which has poles at the eigenvalues: $$g_N(z) ={1 \over N} \sum_{i=1}^N {1 \over z – \lambda_i}.$$

We can relate these two encodings by thinking of our indicator function as the expectation of $1/(z – \lambda)$ relative to the probability distribution we defined: $$ g_N(z) = \int_{-\infty}^\infty { \rho_N(\lambda) \over z – \lambda} d\lambda.$$

It’s then natural to take the $N \to \infty$ limit, and arrive at the Stieltjes transform of our density, $$ g(z)= \int_{-\infty}^\infty { \rho(\lambda) \over z – \lambda} d\lambda,$$ where our discrete distribution $\rho_N$ has now been replaced by the density, $\rho$, that we’re after.

Two questions arise at this point:

1. The function $g(z)$ encodes the density, but we don’t know the density! So how do we compute $g$ without that knowledge?
2. If we somehow manage to compute $g(z)$ by some other means, how do we extraction the density from it?

To answer the first, we return to the finite dimensional case and notice that $g_N(z)$ is the trace of the resolvent:

$$g_N(z) = {1 \over N} \sum_{i=1}^N {1 \over z – \lambda_i} = {1 \over N} \tr((zI – X_N)^{-1}).$$ So, by analogy, $$g(z) \propto \tr((zI – X)^{-1}).$$ The general strategy will then be to compute the finite dimensional trace, and extend somehow to the infinite-dimensional case. The main point is that this doesn’t require direct knowledge of the density.

So now for our second question: Suppose we have $g(z)$, how do we extract the density?

This is where complex analysis comes in. Notice that our function $g(z)$ is undefined on the real line in the ordinary sense, due to the $(z – \lambda)$ in the denominator. However, it’s well-behaved (analytic i.e. complex differentiable – we haven’t shown this) everywhere else, in particular above and below the line. The real line in this case is called a “jump”, because the value of the function “jumps” across it.

We can’t evaluate our function on the real line, but we can consider taking limit values of our function from below the real line: $$g(x – 0 i) \triangleq \lim_{\veps \to 0^+} g(x – \veps i),$$ where by $0^+$ I mean taking the limit from the right.

The Sokhotski–Plemelj theorem then relates this limit to the value of our density on the line. In particular, it says that, for any interval $[a,b]$ containing 0, $$ \im\left[ \lim_{\veps \to 0^+} \int_{a}^b {f(x) \over x – \veps i} dx \right] = \pi f(0).$$

To apply this to our problem, we first extend the integration interval to $[-\infty, \infty]$. We then have $$ g(x – 0 i) = \lim_{\veps \to 0^+} \int_{-\infty}^\infty {\rho(t) \over t – \veps i – \lambda} dt = \lim_{\veps \to 0^+} \int_{-\infty}^\infty {\rho(t – \lambda) \over t – \veps i} dt.$$

This is just the form needed to apply the theorem, so we get $$ \im{g(x – 0 i)} = \pi \rho(x).$$

So, the recipe is:

1. Use the resolvent to determine $g(z)$.
2. Evaluate directional limits of $g(z)$ on the real line.
3. Use the S-P theorem to read off the density from those limits.

Miscellaneous Statements

Below are some elaborations on miscellaneous statements from the book chapter.

They state the definition of the finite dimensional transform $$g_N(z) = \int_{-\infty}^\infty {\rho_N(\lambda) \over z – \lambda} d\lambda,$$ then say that’s well behaved at infinity, and give a series expansion. The expansion can be derived by using the series representation $$ {1 \over z- \lambda} = {1 \over z}{1 \over 1 – {\lambda \over z}} = {1 \over z} \sum_{k=0}^\infty {\lambda^k \over z^k} = \sum_{k=0}^\infty {\lambda^k \over z^{k+1}}.$$ Then $$g_N(z) = \int_{-\infty}^\infty \rho_N(\lambda)\sum_{k=0}^\infty {\lambda^k \over z^{k+1}} d\lambda = \sum_{k=0}^\infty {1 \over z^{k+1}} \int \rho_N(\lambda) \lambda^k d\lambda = \sum_{k=0}^\infty {1 \over z^{k+1}} {1 \over N} \tr(A^k).$$

The importance of infinity is that $z$ will then be beyond all the eigenvalues, so there is no risk of $g_N(z)$ blowing up.

Later, they mention that “knowledge of $g(z)$ near infinity is equivalent to knowledge of all the moments of $A$. What they mean by this is in the moment expansion above.

$$\blacksquare$$