Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!


Ok, here is a very basic example:

n = 1000;
m = 1001;
μ = 0;
σ = 1/Sqrt[n];
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.

share|improve this question
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
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 := – R. M. Mar 21 '14 at 15:54

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.