How to construct the POE ensemble in julia

Question

I'm having a trouble in building the POE ensemble in julia. I am following this paper and part of this other paper.

In julia, I calculate:

X = randn(dim, dim)
Q, R = qr(X)
Q = Q*diagm(sign(diag(R)))
ij = (irealiz-1)*dim
phases_ens[1+ij:ij+dim] = angle(eigvals(Q))

where dim is the matrix dimension and irealiz is just and index for the total number of realizations.

I am interested in the phases of Q, since I want that Q be an orthogonal matrix with the appropriate Haar measure. If dim=50 and the total number of realization is 100000, and since I am correcting Q, I should expect a flat phases_ens distribution. However, I obtain a flat distribution except a peak at zero and at pi. Is there something wrong with the code?

Answer 1

The code is actually correct, you just have the wrong field

The eigenvalue result is true for unitary matrices (complex entries); based on the code from section 4.6 of the Edelman and Rao paper, if you replace the first line by

X = randn(dim, dim) + im*randn(dim, dim)

you get the result you want.

Orthogonal matrices (real entries) behave slightly differently (see remark 1, in section 3 of this paper):

• when dims is odd, one eigenvalue will be +1 or -1 (each with probability 1/2), all others will occur as conjugate pairs.
• when dims is even, both +1 and -1 will be eigenvalues with probability 1/2, otherwise there are no real eigenvalues.

(Thanks for the links by the way: I wasn't aware of the Stewart paper)