I encountered a weird problem when trying to study a few properties of eigenvectors of matrices sampled from the Gaussian symplectic ensemble (GSE). I have encountered this while trying to understand the connection between random matrices and quantum chaos via entanglement entropy, so the motivation is from Physics.
Problem: We want to find the entanglement entropy of the eigenstates of GSE matrices. Given a vector $|\psi \rangle$, we first calculate the reduced density matrix and then calculate entanglement entropy using its eigenvalues. The relevant Mathematica code is (this code finds the reduced density matrix and has been tried and tested many times before)
rdm[\[Psi]ABC_?VectorQ, {dA_Integer /; dA >= 1, dB_Integer /; dB >= 1,
dC_Integer /; dC >= 1}] /; Length[\[Psi]ABC] == dA*dB*dC :=
With[{P = ArrayReshape[\[Psi]ABC, {dA, dB, dC}]},
Flatten[Transpose[P, {1, 3, 2}] .
ConjugateTranspose[P], {{1, 2}, {3, 4}}]]
traceout[\[Psi]_?VectorQ, d_Integer /; d >= 1] /;
Divisible[Length[\[Psi]], d] := rdm[\[Psi], {1, d, Length[\[Psi]]/d}]
traceout[\[Psi]_?VectorQ, d_Integer /; d <= -1] /;
Divisible[Length[\[Psi]], -d] :=
rdm[\[Psi], {Length[\[Psi]]/(-d), -d, 1}]
As there is randomness, we average over many realisations. One can simply do it using the following code (it calculates entanglement entropy of all the eigenstates of a randomly generate GSE matrices and iterate it over savg number of realization)
eestategse = ParallelTable[
{eval, eigveclist} =
Eigensystem[
RandomVariate[
GaussianSymplecticMatrixDistribution[2^(\[CapitalNu] - 1)]]];
Table[
\[Rho]red = traceout[eigveclist[[i]], 2^(\[CapitalNu]/2)];
eigval = DeleteCases[Chop[Eigenvalues[\[Rho]red], 10^-7], 0];
-Total[eigval Log[eigval]], {i, 1, Length[eigveclist]}], {k, 1,
savg}];
The exact result for the eigenstate entanglement for GOE, GUE and GSE type matrices are available in the paper (see Eq. 47,48 and 53). They have a table showing results for some small-dimensional matrices, and one can easily verify it using the above code.
I could verify the numerical results against the analytical results, and they are in agreement for GOE and GUE-type random matrices. However, for the GSE case, I found that the numerical and analytical results do not agree, and the numerical results instead agree with the GUE results. I have done almost everything, like increasing more number of realisations and checking explicitly if the matrices generated are indeed correct or not, but the results are robust.
Based on my understanding, my hunch is that it might have to do with the fact that we use $2 N \times 2N$ dimensional complex matrix representation of $N \times N$ quaternionic symplectic matrices (GSE), hence eigenvectors are kind of behaving like that of complex random matrices (GUE). Also, since the eigenvalues are degenerate so the two eigenvectors that we obtain might both contribute "via some linear combination" to the entanglement entropy, which I do not know (or there is a way to obtain the proper eigenvector using both of them). But I am sure that indeed they are a correct GSE matrix, as one can find the eigenvalue spacing distribution of the same, and it matches with the Wigner-Dyson for GSE, so no problem with the matrices. From the Physics side, one can also check that it is invariant under $\mathcal{T}$ such that $\mathcal{T} H_{GSE} \mathcal{T}^{-1} = H_{GSE}$ where $H_{GSE}$ is a GSE matrix and $\mathcal{T}$ is an antiunitary operator such that $\mathcal{T}^{2}= -1$. I have checked all these properties, and they hold.
I would thus appreciate any insight on this that can help me \emph{restore} the properties of GSE eigenvectors, because for now it seems that what Mathematica or, in fact, any other programming language gives are not behaving like a quaternionic vector. In case I am missing something crucial here regarding these matrices, please let me know as well. To my end, I also tried to look in the literature regarding any such numerical calculation, but could not really find any and hence the question (as of now).
