# What is the fastest way to obtain the eigenvalues of a Wishart matrix?

I would like a fast method of creating a sample of random numbers which corresponds to the eigenvalues of a Wishart matrix:

For M>N the eigenvalues \lambda_i are given by the jpd

Where $K_N$ is a normalisation constant and $\beta = 1,2$ depends on whether $X$ is real or complex. Does anyone have any suggestions for a fast way of implementing this? My goal is to play with multiple realisations of eigenvalues representing 100x100 matrices.

EDIT:

Ok, here is a very basic example:

n = 1000;
m = 1001;
μ = 0;
σ = 1/Sqrt[n];
AbsoluteTiming[
a = RandomVariate[NormalDistribution[μ, σ], {n, m}];
wish = a.a\[ConjugateTranspose];
ewish = Eigenvalues@wish;
]
ewishPlot = Histogram[ewish, {0.1}, "PDF"];

ηp = n σ^2 (1 + Sqrt[η])^2;
ηm = n σ^2 (1 - Sqrt[η])^2;
η = m/n;
ρmp[λ_] := 1/(2 π n σ^2 λ)Sqrt[(ηp - λ) (λ - ηm)];
ρmpPlot = Plot[ρmp[λ], {λ, 0, 4}];
Show[ρmpPlot, ewishPlot]


I guess speed won't be a problem but I would still very much like to know the best way of doing this.

-
Do you have a slow method now? Wishart matrices are trivial to construct and Eigenvalues on 100x100 machine number matrices should be very fast as well... – R. M. Mar 21 '14 at 14:04
Have a look here. I would've done that probably differently now, but that could be a start. – Leonid Shifrin Mar 21 '14 at 14:31
In fact, about 6 years ago I was generating a bunch of these for betas 2 and 4 (IIRC), to study eigenvalue statistics in the context of applications to QCD (Dirac operator in Random Matrix Theory), so I could probably dig up some of that code, but what is in the book can also be a starting point. – Leonid Shifrin Mar 21 '14 at 14:34
Thanks very much. This looks great! – user12876 Mar 21 '14 at 15:22
Because it is implemented more generally than what you've done here. You can use a non-diagonal Sigma and it uses MultinormalDistribution to generate the data. In any case, the core part of it is implemented the same way as yours, so I don't see a reason to look for a different way. It's probably as efficient as it can be. – R. M. Mar 21 '14 at 18:11

Since all the entries of $X$ are i.i.d. normal, it can be converted into a band diagonal matrix with $\chi$ distributed components via elementary reflections. A paper discusses about this is "Matrix Models for Beta Ensembles" and can be found in the following link: