Is tridiagonal reduction the current best practice to compute eigenvalues of random matrices from the Gaussian ensembles (GOE, GUE, GSE)?

Question

I have tried to compute the eigenvalues of random matrices of the GOE ensemble, using MATLAB.

Such matrices of size $n * n $ can be obtained easily, symmetrizing matrices whose elements follow the standard normal distribution: $H = (M + M^T )/ 2 $.

In MATLAB code for instance:

H = randn(n); 
M = (H+H')/2;
eigenvalues = eig(M);

This naive approach is the one I found in most documents. However, it looks like a more efficient method was developed some years ago by A. Edelman and others, which makes use of a matrix factorization (Householder factorization) to obtain a tridiagonal matrix whose eigenvalues are following the same distribution than the GOE. A description of this technique can be found here: http://www-math.mit.edu/~edelman/publications/random_matrix.pdf

The form of this tridiagonal matrix is the next one:

\begin{bmatrix} G_n & \chi_{n-1} & & & & \\ \chi_{n-1} & G_{n-1} & \chi_{n-2} & & & & \\ & \chi_{n-2} & G_{n-2} & & & \\ & & & \ddots & \ddots & \chi_{1} \\ & & & & \chi_{1} & G_1 \\ \end{bmatrix}

where G follow a normal distribution and $\chi_i$ follow chi square distributions. Such a matrix in MATLAB can be obtained with the next lines:

a = sqrt(chi2rnd([n:-1:1]))'; 
H = spdiags(a, 1, n, n) + spdiags(randn(n,1)/sqrt(2), 0, n, n);
M = (H+H')/sqrt(2);
eigenvalues = eig(M);

Eigenvalues from such a matrix should be easier to get than the ones of the initial problem (the author comes up with $O(n^2)$ time vs. $O(n^3)$ time). When I plot some histograms of the eigenvalues obtained with the two methods, I get results in good agreement.

Shall this method be considered the current best practice to compute the eigenvalues of matrices from the GOE ? What do people do in practice ?

Answer 1

Q: What do people do in practice ?

The "practice" will be field specific, but in most of the physics applications I am aware of one needs not only the eigenvalues but also the eigenvectors. For example, the GOE may be used to model the Hamiltonian of a small metal grain or a semiconductor quantum dot, and one wants to calculate transport properties such as the conductance of that object. The eigenvalues of the Hamiltonian do not give sufficient information.

So at least in that context one will want to generate the full matrix ensemble, not only the eigenvalue distribution.

Answer 2

Best practice for what please?: speed, accuracy, numerical stability, memory load, scaling on hardware parallel architectures...? Do you need all eigenvalues or only some of them? Do you need their multiplicities? Do you also need the eigenvectors as point out by Carlo?

I could not find a single eigenvalue/eigenvector algorithm tailored to Gaussian matrices. It seems the Gaussian structure is quite useless from the numerical algebra point of view.

So, we are left with the standard eigenproblems for symmetric matrices for the GOE. Reduction to tridiagonal form is common practice for such problems, this is probably what Matlab is doing (and no doubt Matlab is the best software for numerical linear algebra).

But 1) there are many different methods to reduce to the tridiagonal form (e.g. we can use Lanczos instead of Householder) 2) then there are many differents methods to get the eigenvalues (e.g. in $O(n\log(n)$) instead of $O(n^2)$ by divide-and-conquer methods) 3) reduction to tridiagonal form is only one approach among many others (e.g. reduction to Hessenberg form then QR, Jacobi........).

Numerical linear algebra is a jungle.