We are interested in calculating the spectral moments of random matrices sampled from the Wishart distribution. Let $N,T$ be positive integers with $T> N$ . Then the quantities in question are defined as $m_p:= 1/N \cdot Tr \left[ c^p \right]$ where $p \in {\mathbb Z}$ and where $c:= 1/T \cdot {\bf X} \cdot {\bf X}^T$ is the sample correlation matrix and ${\bf X}:= \left( X_{i,t} \right)_{i=1,t=1}^{N,T}$ is a random matrix with iid entries being sampled from a normal distribution with mean zero and variance one.
Let us define ${\mathfrak P}_{N,T} = \frac{\pi ^{\frac{N^2}{2}-\frac{1}{2} (N-1) N} 2^{-\frac{N T}{2}}}{\left(\prod _{j=0}^{N-1} \Gamma \left(\frac{N-j}{2}\right)\right) \prod _{j=0}^{N-1} \Gamma \left(\frac{T-j}{2}\right)}$ and ${\mathfrak a} := (T-N-1)/2 $. Then, as it turns out, those spectral moments adopt a following neat closed form expression:
\begin{eqnarray} m_p = \frac{{\mathfrak P}_{N,T}}{N \cdot T^p} \cdot 2^{\binom{N+1}{2} +N {\mathfrak a} +p} \cdot \sum\limits_{J=1}^N \left. \prod\limits_{\xi=1}^N ( \partial_{{\bar A}_\xi} + \cdots + \partial_{{\bar A}_N} )^{{\mathfrak a} + \delta_{\xi,J} p} \cdot \prod\limits_{1 \le \xi < \eta \le N} ( \partial_{{\bar A}_\xi} + \cdots + \partial_{{\bar A}_{\eta-1}} ) \cdot \prod\limits_{\eta=1}^N \frac{1}{{\bar A}_\eta} \right|_{{\bar A}_\eta=\eta} \tag{1} \end{eqnarray}
Equation $(1)$ follows from the knowledge of the joint distribution of eigenvalues in the Wishart ensemble (see
Harnad, John (ed.), Random matrices, random processes and integrable systems, CRM Series in Mathematical Physics. Berlin: Springer (ISBN 978-1-4419-9513-1/hbk; 978-1-4419-9514-8/ebook). xviii, 524 p. (2011). ZBL1215.15002-Theorem 7.1.2 "Wishart Distribution" section 7.1.1 page 445).
That distribution involves and integral over the orthogonal group -- the so called Itzykson-Zuber-Harish-Chandra integral -- which is unsolved. However, in the case when the underlying correlation matrix in the ensemble in question is an identity matrix the aforementioned integral reduces to a unity and a neat closed form expression is obtained, an expression which can be then interpreted as the distribution of energies in the so called Dyson Brownian motion model (see
Dyson, F. J., A Brownian-motion model for the eigenvalues of a random matrix, J. Math. Phys. 3, 1191-1198 (1962). ZBL0111.32703.).
The evaluation of the spectral moments reduces then to evaluating a multi-dimensional integral over the Weyl chamber, a task that can be then translated to evaluating an action of a multivariate differential operator on a product of reciprocals as explained in here.
After this lengthy introduction we tried to put our hands on equation $(1)$ and see if we can evaluate it in closed form. We took $N=2$, we expanded the multivariate differential operator and carried out its action analytically and we arrived at a following expression:
\begin{eqnarray} &&m_p = \frac{{\mathfrak P}_{2,T}}{2 \cdot T^p} \cdot 2^{3+2 {\mathfrak a}+p} \cdot ({\mathfrak a})! ({\mathfrak a}+p)! \cdot \left[ \right. \\ && % \left. \sum\limits_{l=0}^{{\mathfrak a}+p} \frac{(1+{\mathfrak a}+p-l)^{({\mathfrak a})}}{1^{({\mathfrak a})}} \cdot \frac{(l+1)}{1^{l+2}} \cdot \frac{1}{2^{2 {\mathfrak a}+p-l+1}} + \right. \\ && % \left. \sum\limits_{l=0}^{{\mathfrak a}} \frac{(1+{\mathfrak a}-l)^{({\mathfrak a}+p)}}{1^{({\mathfrak a}+p)}} \cdot \frac{(l+1)}{1^{l+2}} \cdot \frac{1}{2^{2 {\mathfrak a}+p-l+1}} \right. \\ && \left. \right] \tag{2} \end{eqnarray}
Then we simply evaluated expression $(2)$ for $p=0,-1,-2$ and for $T=2+1-2 p- 2 \xi$ where $\xi=0,\cdots, 3$ and we discovered the following pattern:
\begin{eqnarray} && m_p = \left\{ \begin{array}{lll} 1 & \mbox{if $p=0$} \\ \frac{T}{T-3} & \mbox{if $p=-1$} \\ \frac{T^2 (-1+T)}{(-5+T)(-3+T)(-2+T)} & \mbox{if $p=-2$} \\ \vdots \end{array} \right. \tag{2a} \end{eqnarray}
The code snippet below verifies the expressions in $(2a)$. We have:
In[65]:= mGamma[NN_, x_] :=
Pi^(NN (NN - 1)/4) Product[Gamma[(2 x - j)/2], {j, 0, NN - 1}];
Pfct[NN_, T_] :=
Pi^(NN^2/2) 2^(-NN T/2)/(mGamma[NN, NN/2] mGamma[NN, T/2]);
ll1 = Table[
With[{T = 2 + 1 - 2 p + 2 xi},
With[{a = (T - 2 - 1)/2},
Pfct[2, T]/(2 T^p) 2^(3 + 2 a + p) (a + p)! a! (
Sum[Pochhammer[1 + a + p - l, a]/Pochhammer[1, a] (l + 1) 1/
2^(2 a + p - l + 1), {l, 0, a + p}] +
Sum[Pochhammer[1 + a - l, a + p]/
Pochhammer[1, a + p] (l + 1) 1/2^(2 a + p - l + 1), {l, 0,
a + 0}])]], {p, 0, -2, -1}, {xi, 0, 3}]
ll2 = Table[
Which[p == 0, 1, p == -1,
With[{T = 2 + 1 - 2 p + 2 xi}, T/(T - 2 - 1)], p == -2,
With[{T =
2 + 1 - 2 p + 2 xi}, ((-1 + T) T^2)/((-5 + T) (-3 + T) (-2 +
T))], True, NaN], {p, 0, -2, -1}, {xi, 0, 3}];
ll1 - ll2
Out[67]= {{1, 1, 1, 1}, {5/2, 7/4, 3/2, 11/8}, {147/20, 27/7, 605/216,
507/220}}
Out[69]= {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}
Now I have two questions. Firstly, can we find a closed form expression for the spectral moment in $(2)$ for arbitrary values of $p = -3,-4,\cdots$ ?
Secondly, can we find similar expressions for the moments for $N > 2 $?
