[proofplan]
We diagonalize the real symmetric matrix $A$ by an orthogonal change of basis. Since $z$ is non-real and all eigenvalues of $A$ are real, each scalar $\lambda_i-z$ is non-zero, so the resolvent $(A-zI_n)^{-1}$ exists and is diagonalized by the same orthogonal matrix. Taking the trace reduces the resolvent trace to the average of the scalar reciprocal terms $(\lambda_i-z)^{-1}$, which is exactly the [Stieltjes transform](/page/Stieltjes%20Transform) of the [empirical spectral distribution](/page/Empirical%20Spectral%20Distribution).
[/proofplan]
[step:Diagonalize $A$ and invert the shifted diagonal matrix]
The matrix $A$ represents a self-adjoint endomorphism of the finite-dimensional real [inner product space](/page/Inner%20Product%20Space) $\mathbb{R}^n$ because $A=A^\top$. By the [Spectral Theorem For Self-Adjoint Maps](/theorems/440), there exist an orthogonal matrix $Q \in \mathbb{R}^{n \times n}$ and a diagonal matrix $\Lambda \in \mathbb{R}^{n \times n}$ such that
\begin{align*}
A = Q\Lambda Q^\top.
\end{align*}
The diagonal matrix is
\begin{align*}
\Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_n).
\end{align*}
Fix $z \in \mathbb{C}\setminus\mathbb{R}$. Since each $\lambda_i \in \mathbb{R}$, we have $\lambda_i-z \neq 0$ for every $i \in \{1,\dots,n\}$. Hence $\Lambda-zI_n$ is invertible over $\mathbb{C}^{n \times n}$, with
\begin{align*}
(\Lambda-zI_n)^{-1}
=
\operatorname{diag}\left((\lambda_1-z)^{-1},\dots,(\lambda_n-z)^{-1}\right).
\end{align*}
Since $I_n = QI_nQ^\top$, we have
\begin{align*}
A-zI_n = Q\Lambda Q^\top - zQI_nQ^\top = Q(\Lambda-zI_n)Q^\top.
\end{align*}
Therefore
\begin{align*}
(A-zI_n)^{-1}
=
Q(\Lambda-zI_n)^{-1}Q^\top.
\end{align*}
[guided]
The point of diagonalizing $A$ is that the inverse of $A-zI_n$ becomes the inverse of a diagonal matrix after changing basis. The hypothesis $A=A^\top$ says exactly that $A$ represents a self-adjoint endomorphism of the finite-dimensional real [inner product](/page/Inner%20Product) space $\mathbb{R}^n$. Therefore the [Spectral Theorem For Self-Adjoint Maps](/theorems/440) applies and gives an orthogonal matrix $Q \in \mathbb{R}^{n \times n}$ and a diagonal matrix $\Lambda \in \mathbb{R}^{n \times n}$ satisfying
\begin{align*}
A = Q\Lambda Q^\top.
\end{align*}
The diagonal matrix is
\begin{align*}
\Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_n),
\end{align*}
where $\lambda_1,\dots,\lambda_n$ are the real eigenvalues of $A$, counted with algebraic multiplicity.
Now fix $z \in \mathbb{C}\setminus\mathbb{R}$. The non-real hypothesis on $z$ is used exactly here: since every $\lambda_i$ is real, no number $\lambda_i-z$ can vanish. Thus the shifted diagonal matrix $\Lambda-zI_n$ is invertible, and its inverse is obtained by inverting each diagonal entry:
\begin{align*}
(\Lambda-zI_n)^{-1}
=
\operatorname{diag}\left((\lambda_1-z)^{-1},\dots,(\lambda_n-z)^{-1}\right).
\end{align*}
Because $Q$ is orthogonal, $Q^\top Q = I_n$ and $QI_nQ^\top = QQ^\top = I_n$. Using $I_n = QI_nQ^\top$, we rewrite the shifted matrix as
\begin{align*}
A-zI_n = Q\Lambda Q^\top - zI_n = Q\Lambda Q^\top - zQI_nQ^\top = Q(\Lambda-zI_n)Q^\top.
\end{align*}
Multiplying $Q(\Lambda-zI_n)^{-1}Q^\top$ on either side by $A-zI_n$ gives the identity matrix, using $Q^\top Q=I_n$. Therefore
\begin{align*}
(A-zI_n)^{-1}
=
Q(\Lambda-zI_n)^{-1}Q^\top.
\end{align*}
[/guided]
[/step]
[step:Compute the resolvent trace from the diagonal entries]
We regard the real matrices $Q$ and $Q^\top$ as complex matrices by scalar extension. Define $D \in \mathbb{C}^{n \times n}$ by $D=(\Lambda-zI_n)^{-1}$, and write $Q_{ij}$ for the $(i,j)$ entry of $Q$. For each $i \in \{1,\dots,n\}$, the $(i,i)$ entry of $QDQ^\top$ is
\begin{align*}
(QDQ^\top)_{ii}
=
\sum_{k=1}^{n}\sum_{\ell=1}^{n} Q_{ik}D_{k\ell}(Q^\top)_{\ell i}.
\end{align*}
Since $D$ is diagonal, $D_{k\ell}=0$ when $k\neq \ell$, so this becomes
\begin{align*}
(QDQ^\top)_{ii}
=
\sum_{k=1}^{n} Q_{ik}D_{kk}Q_{ik}.
\end{align*}
Summing over $i$ and interchanging the two finite sums gives
\begin{align*}
\operatorname{tr}\left((A-zI_n)^{-1}\right)
=
\operatorname{tr}(QDQ^\top)
=
\sum_{k=1}^{n}D_{kk}\sum_{i=1}^{n}Q_{ik}^2.
\end{align*}
The columns of the orthogonal matrix $Q$ have Euclidean norm $1$, so $\sum_{i=1}^{n}Q_{ik}^2=1$ for every $k$. Hence
\begin{align*}
\operatorname{tr}\left((A-zI_n)^{-1}\right)
=
\sum_{k=1}^{n}D_{kk}
=
\operatorname{tr}(D)
=
\operatorname{tr}\left((\Lambda-zI_n)^{-1}\right).
\end{align*}
Since $(\Lambda-zI_n)^{-1}$ is diagonal, its trace is the sum of its diagonal entries:
\begin{align*}
\operatorname{tr}\left((A-zI_n)^{-1}\right)
=
\sum_{i=1}^{n}\frac{1}{\lambda_i-z}.
\end{align*}
[/step]
[step:Identify the scalar sum with the Stieltjes transform]
By the definition of the [empirical spectral distribution](/page/Empirical%20Spectral%20Distribution),
\begin{align*}
\mu_A = \frac{1}{n}\sum_{i=1}^{n}\delta_{\lambda_i}.
\end{align*}
Let $f_z: \mathbb{R}\to \mathbb{C}$ be the function defined by $f_z(x)=(x-z)^{-1}$ for every $x\in\mathbb{R}$. By the definition of the [Stieltjes transform](/page/Stieltjes%20Transform) and the defining property of the Dirac measure, we have
\begin{align*}
m_{\mu_A}(z)
&= \int_{\mathbb{R}} f_z(x)\, d\mu_A(x),
\end{align*}
and hence
\begin{align*}
m_{\mu_A}(z)
&= \frac{1}{n}\sum_{i=1}^{n}\int_{\mathbb{R}} f_z(x)\, d\delta_{\lambda_i}(x),
\end{align*}
so
\begin{align*}
m_{\mu_A}(z)
&= \frac{1}{n}\sum_{i=1}^{n}\frac{1}{\lambda_i-z}.
\end{align*}
Combining this identity with the trace computation gives
\begin{align*}
m_{\mu_A}(z)
=
\frac{1}{n}\operatorname{tr}\left((A-zI_n)^{-1}\right).
\end{align*}
Since $z \in \mathbb{C}\setminus\mathbb{R}$ was arbitrary, the formula holds for every non-real spectral parameter $z$.
[/step]
