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

Question

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

definition of Wishart matrix

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

eigenvalue distribution

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
What you have looks pretty good to me (this is how I've done it). The alternate that I know of is "cleaner", but very slow — Needs["MultivariateStatistics`"]; eigs = Eigenvalues@RandomReal[WishartDistribution[σ IdentityMatrix@n, m]] You might want to define your functions with :=... see here for the difference between = and := — rm -rf, Mar 21 '14 at 15:54
Why is that slower? I would have expected that if there is a built in function called "WishartDistribution" it would be better to use it than build one from scratch. — user12876, Mar 21 '14 at 16:09
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

score 3 · Answer 1 · edited May 10 '20 at 15:16

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 this arXiv preprint.

The Wishart matrix constructed in this way is tridiagonal and can be represented as SparseArray in Mathematica. The eigenvalues can be obtained efficiently by using Eigenvalues with method "Banded".

score 2 · Answer 2 · answered Oct 03 '18 at 04:33

One can now use MatrixPropertyDistribution[] along with WishartMatrixDistribution[] to generate the eigenvalues of a random Wishart matrix. For instance,

RandomVariate[MatrixPropertyDistribution[Eigenvalues[X], 
              X \[Distributed] WishartMatrixDistribution[4, IdentityMatrix[4]]], 5]

will generate $5$ sets of eigenvalues of random $4\times4$ Wishart matrices.

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

2 Answers2