Symmetry lost in calculating eigenvectors with julia

Question

Can someone point out why julia loses symmetries when diagonalizing certain types of matrices? (In the following I will ignore normalization constants.) I've been trying to solve the following Floquet matrix:

U[i,i] = exp(-i^2im/N),

with N the dimension of the matrix and other components are zero. Clearly, this is the "time-evolution" of the Hamiltonian

H = p^2/2.

This H is symmetric under parity and so is U in the site basis Us = Udft'U*Udft (Udft is the discrete Fourier matrix s.t. Udft[m,mp] = sqrt(N)^-1 exp(i*j*1im/N) see below), that is, one can check that

Jp*Us - Us*Jp = 0,

where

Jp[i,j] = \delta_{i,N-j+1}

is the space inversion matrix. The eigenstates, however, do not fulfill parity. If vs is an eigenstate of Us, then

Jp*vs = \pm vs,

which does not happen for the numerical results given by julia. It's a little bit bizarre because for low dimensionality, say N=11 there is no problem, but if I go to, say N=1001 then trouble begins to apppear. (For certain conditions, I want N to be odd. The reason is that I have the particle constrained to move over the unitary circle and I want that the sites are symmetric around the angle zero.) To diagonalize I use Julia 1.2.0 and

 LinearAlgebra.eigen(Us).

Addendum: Thanks for SGJ for pointing out at apparent error in the DFT. To construct the matrix I do

M = div(N,2)
m = 1
for ii in -M:M
    mp = 1
    for jj in -M:M
        U[m,mp] = exp(2.0im*pi*ii*jj/Nsites) #                             
        mp += 1
    end
    m += 1
end

so that positive and negative frequencies enter in the DFT and are chosen symmetric around momentum zero from -N/2 to N/2 where integer division is used.

Answer 1

First, roundoff errors will generally break symmetries slightly, so small discrepancies from mirror symmetry (relative errors norm(Jp*vs - ±vs)/norm(vs) < 1e-12) should not be a concern.

Second, the p^2 operator's energy eigenvalues are doubly degenerate (left- and right-going waves have the same energy), and parity of degenerate states is not guaranteed by symmetry. Any linear combination of the degenerate even and odd eigenfunctions (= standing waves) is also an eigenfunction, so you can make non-symmetric eigenfunctions (e.g. a left-going wave rather than a standing wave). When you have a degeneracy, the correct statement is that you can choose an eigenfunction basis that is also a parity eigenfunction. However, Julia's eigen operation (the standard LAPACK algorithm) basically chooses a "random" basis for a degenerate eigenspace, so it will not necessarily be the one you want.

Third, it doesn't look like you're using the discrete Fourier transform (DFT) correctly. You're essentially using the DFT to express the second-derivative operator (p^2) in a planewave (momentum-space) basis, but your posted answer is forgetting about aliasing — you really need to have both positive and negative frequencies present. That is, the frequency of half the terms should be proportional to i-1 (due to Julia's 1-based indexing) and the other half proportional to N+1-i. This is explained in detail in many sources, e.g. these notes.