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.

Many thanks in advance!

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... – rm -rf 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. – rm -rf 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:
will generate $$5$$ sets of eigenvalues of random $$4\times4$$ Wishart matrices.