# Why is generating normalized random $1000 \times 1000$ matrices and plotting the eigenvalues so slow?

For each of the distributions $$N(0,1)$$ and $$\pm 1$$ equal probability and for each of $$N \in \{5,10,20,50,100,200,1000\},$$ I want to generate an $$N \times N$$ matrix with entries chosen from the distribution, normalize the matrix by dividing all entries by $$\sqrt{N},$$ then plot the eigenvalues. Finally, repeat the whole process but with random symmetric matrices. The eigenvalues will be real, so in this case the plot becomes a density histogram.

The general and symmetric case already take about 10 seconds to run without $$N = 1000,$$ but when I add $$1000$$ it hasn't finished running after several minutes. Even $$N = 500$$ took too long. What's going on? Code is shown below:

cp[A_] := ComplexListPlot[Eigenvalues[A]]
rp[A_] := Histogram[Eigenvalues[A], {1/Sqrt[Length[A]]}, "PDF"]
randSym[dist_, n_] := Module[{diag, up}, diag = RandomVariate[dist, n];
up = RandomVariate[dist, Binomial[n, 2]];
StatisticsLibraryVectorToSymmetricMatrix[up, diag, n]]

rMat[dist_, n_] := RandomVariate[dist, {n, n}]
sizes = {5, 10, 20, 100, 200}
dists = {TransformedDistribution[-1 + 2 x,
x \[Distributed] BernoulliDistribution[1/2]],
NormalDistribution[0, 1]}
Outer[gPlot, dists, sizes]
Outer[symPlot, dists, sizes]


Plots look as expected for $$N \le 200,$$ so runtime is the only issue. Plotting, scaling, and entry filling are linear with the amount of data, so I suspect the eigenvalue computation needs to be optimized. According to this answer, it should only take 5 seconds for $$N = 1000,$$ so where did I go wrong?

• Do you need exact results? NormalDistribution produces machine-precision numbers (by default), but other distributions you use don't (returning exact values), and this results Eigenvalues computing massively complicated roots. If numerics are sufficient, you can accomplish this, for example, by replacing Eigenvalues[A] with Eigenvalues[N[A]]. This speeds computation immensely. Feb 23, 2023 at 5:11

You are probably unintentionally computing exact eigenvalues since your RandomVariates from TransformedDistribution are exact numbers in this case. The toy example below corresponds to the problem you're having:

With[
{n = 40,
dist = EmpiricalDistribution[{-1, 1}]},
Eigenvalues[
StatisticsLibraryVectorToSymmetricMatrix[
RandomVariate[dist, Binomial[n, 2]],
RandomVariate[dist, n],
n]/Sqrt[n]]] // First

(* Root[-150375024359015916240896 + 1441199744254536367210496 # +
2325218250072738087567360 #^2 - 13017297010673087538003968 #^3 +
564877696112481338916864 #^4 + 29671492037768896734822400 #^5 -
17802351003049187799465984 #^6 - 18093175893422395405369344 #^7 +
15880355366435254528114688 #^8 + 4785459350035091689570304 #^9 -
6289499837890258177884160 #^10 - 536516394575471963013120 #^11 +
1449836412440797431988224 #^12 - 14137372818718179983360 #^13 -
217935992237474696921088 #^14 + 12515845805302628745216 #^15 +
22753565305512700936192 #^16 - 1855539065363101646848 #^17 -
1715940016012178489344 #^18 + 158028501898907090944 #^19 +
95858785732624121856 #^20 - 9016655135417827328 #^21 -
4029693175314907136 #^22 + 364081625705611264 #^23 +
128561360111075328 #^24 - 10652592683483136 #^25 -
3119691379490816 #^26 + 227617992572928 #^27 + 57350599718912 #^28 -
3538368765952 #^29 - 789882882048 #^30 + 39398755840 #^31 +
7985353088 #^32 - 304575488 #^33 - 57238656 #^34 + 1543744 #^35 +
274056 #^36 - 4584 #^37 - 782 #^38 +
6 #^39 + #^40& , 40, 0]/(2 Sqrt[10]) *)


That's just the first eigenvalue for this matrix in its exact polynomial root form, and handling these roots gets increasingly complicated as the input matrix size grows.

If you are only interested of finite-precision numerical approximations of eigenvalues, probably the easiest way to work around computing exact solutions is to provide input at as machine-precision numbers, which can be accomplished by replacing Eigenvalues[matrix] with Eigenvalues[N[matrix]]. This can provide a dramatic speedup.

• The distributions aren't accepted once I do that: "Nonatomic expression expected at position 2 in \ Outer[gPlot,dists,sizes]" Feb 23, 2023 at 16:27
• @Displayname sounds like you have replaced dists with something that is not a list. Outer expects it to be a list... Feb 23, 2023 at 17:10
• It's still a list. I didn't touch anything, but I came back to rerun all the cells and now there's a new error: "Positive machine-sized integer or Infinity expected at position 3 in \ Outer[gPlot,{TransformedDistribution[-1+2\ [FormalX],[FormalX]\ [Distributed]BernoulliDistribution[1/2]],NormalDistribution[0,1]},\ sizes]." The biggest red flag is that dists is no longer blue, I have no idea why Feb 23, 2023 at 17:17
• Ok, everything is working now. Bizarre for the error to appear and disappear. Feb 23, 2023 at 17:20