# Are Python/MATLAB/Mathematica numerical eigenvectors affected by eigenvalue degeneracies outside region of calculation?

I have a discretized 2D mesh over which I calculate eigenvalues and eigenvectors of some Hermitian 2 x 2 matrix at each point along a closed loop parameterized by parameter t. The eigenvectors are usually complex-valued. I then calculate a quantity z(t) using the eigenvectors. The complex-valued number z(t) is calculated by taking the dot product of the first eigenvector with the derivative of the second eigenvector. The derivative is calculated using a standard central difference scheme along the loop. When z(t) is plotted on the complex plane - i.e. (Re z(t),Im z(t)) - it is a smooth curve.

However, sometimes, there might be a point in the 2D space with an eigenvalue degeneracy. The issue is that the complex loop z(t) ends up a non-smooth, jagged mess even if the eigenvalues are non-degenerate along the loop over the mesh. From the physics problem this code models, this is NOT expected. We expect smooth z(t) as long as there is no eigenvalue degeneracy along the loop. The figure at the end illustrates the context.

Therefore, I am confused. Because, numerically, it's as if the eigenvectors along the loop feel the presence of the eigenvalue degeneracy that occurs elsewhere. It seems to be a numerical artifact, but I don't understand what's going on. I see this behavior in Python/MATLAB/Mathematica implementations of the code. So, does anyone have insight on what might be going on? Is this a universal feature of numerical eigen solvers? What happens at a fundamental level in these solvers that causes this? Or am I wrong in correlating the degeneracy with the non-smooth complex loop? This behavior is observed in different matrix models (even for a pair of eigenvectors from a matrix larger than 2 x 2). Any insight would be appreciated!

EDIT: schematic of calculation

• n_segs is the number of pieces to discretize the loop into (=250 usually)

• kloop is an array of 2D coordinates specifying the path in the (x,y) domain. It is generated by function parameterize_loop() given below

• Hamiltonian() returns the square matrix at given (x,y)

• sort_eigens(), defined below, sorts eigenvectors based on eigenvalue magnitude ordering

• finite difference is calculated in dphim

• wfc_dot(), defined below, takes the dot product of the two vectors (complex conjugate of one vector and the finite difference of the other vector, per standard quantum mechanical inner product procedure)

• cloop is an array of complex numbers that I plot in the figure on the right; i.e. z(t).

 for j in np.arange(1,n_segs):

H = Hamiltonian(kloop[j-1][0],kloop[j-1][1])
eval, psi0 = LA.eig(H)
[eval,psi0]=sort_eigens(eval,psi0);

H = Hamiltonian(kloop[j][0],kloop[j][1])
eval2, psi1 = LA.eig(H)
[eval2,psi1]=sort_eigens(eval2,psi1);

H = Hamiltonian(kloop[j+1][0],kloop[j+1][1])
eval, psi2 = LA.eig(H)
[eval,psi2]=sort_eigens(eval,psi2);

psi0 = psi0[:,band2]
psi2 = psi2[:,band2]
dphim = (psi2-psi0)*n_segs/(4*np.pi);

phip = psi1[:,band1]

cloop[j-1] = wfc_dot(phip,dphim)

cloop[-1]=cloop[0] # this ensures the loop is closed; not really necessary


   def sort_eigens(eigenValues,eigenVectors):
idx = eigenValues.argsort()
vals=eigenValues[idx]
vecs=eigenVectors[:,idx]
for q in range(len(vecs)):
v = vecs[:,q]
non_zero_elems = v[np.nonzero(v)]
v = v/(non_zero_elems[0,0])
return [vals,vecs]


    def wfc_dot(wf1,wf2):
wf1 = np.squeeze(np.asarray(wf1))
wf2 = np.squeeze(np.asarray(wf2))
return np.dot(wf1.flatten().conjugate(),wf2.flatten())


    def parameterize_loop(kx0,ky0,r,n_segs): # parameterizes a counterclockwise loop of radius r centered at (kx0,ky0)
a = np.linspace(0.0, 2*np.pi, num=n_segs+1) # time goes from 0 to 2pi
a = np.append(a,a[1]) # to ensure continuity in subsequent winding loop calculations
k_loop=np.zeros((n_segs+2,2))
for i in range(n_segs+2):
k_loop[i,:]=np.array([kx0 + r*np.cos(a[i]), ky0 + r*np.sin(a[i])])
return k_loop

• Do you have repeated eigenvalues that might split into isolated eigenvalues (or vice versa)? When that happens, there's no guarantee that the corresponding eigenvectors will change in a continuous fashion. Jan 27, 2022 at 17:40
• @BrianBorchers There are no repeated eigenvalues along the path taken. However, this issue arises when a single point elsewhere has the same eigenvalue. Jan 27, 2022 at 17:55
• How are your eigenvalues being computed (e.g. by some specific LAPACK routine or specific numpy routine)? When you use a finite difference approximation for the derivative how do you select the step size for that finite difference? Jan 27, 2022 at 18:42
• How close to having repeated eigenvalues do you get when things blow up? Is the relative difference large or small relative to floating point precision? (single or double precision?) Jan 27, 2022 at 18:44
• In light of my answer, you should probably use LA.eigh(), not eig() Jan 27, 2022 at 19:07

At a glance I too would expect this problem to be pretty well behaved, because Hermitian eigendecomposition itself is pretty well behaved (real eigenvalues, orthogonal eigenvectors). One place where things might go wrong, is calling a routine/algorithm that is unaware of the Hermitian-ness of your input. In LAPACK, the routine you'd probably want is ZHEEV, which does offer guarantees you'd expect on a Hermitian eigendecomposition: orthogonal eigenvectors, real eigenvalues, even sorted for you. You can contrast this with ZGEEV, which can in principle compute the eigendecomposition just as well, but doesn't exploit the structure and doesn't offer the guarantees.

ZHEEV documentation.

ZGEEV documentation.

To clarify, the guarantees are from the problem properties, not the algorithm .. so if you accidentally use ZGEEV instead of ZHEEV, your eigenvectors will still be orthogonal but probably not to the same precision, your eigenvalues will be real but probably with some small/roundoff-level imaginary part, etc. They certainly won't be sorted because ZGEEV has no business trying to apply a sort to eigenvalues that it presumes to be complex. It's also entirely possible that the two algorithms have different sensitivity to degeneracy, because they are just using different transformations under the hood.

The sorting bit seems particularly noteworthy, as I expect you'd need some kind of "tracking"/continuity of the eigenpairs as you vary the t-parameter, otherwise you might inadvertently compute dot(x1,dx2) instead of dot(dx2,x1), or whatnot. Or worse yet, you might inadvertently finite difference x1(t+dt) with x2(t-dt), because they got scrambled/swapped during eigensolution.

So I guess that's what I would dig into first, is making sure that you are informing whatever higher level programming environment you are using, about these symmetry/Hermitian properties so that it can dispatch to the right underlying algorithms. This might take the form of extra arguments ('symm=true', perhaps), or using special Hermitian-only matrix containers/types? Hard to say for certain. And make sure that you are really tracking/selecting the correct eigenvectors when you are differencing/dotting.

• Thank you for the answer. I’ll look into it. It might take me a few days to get back to you. Jan 27, 2022 at 17:53
• Upon rereading that the eigenproblems in question are 2x2, you might even consider solving it using non-LAPACK machinery .. you can eigendecompose such a matrix using just a single Jacobi/double-sided-Givens rotation. Hermitian-ness can be built directly into this calculation and sorting can be performed as a postprocess. Jan 27, 2022 at 18:55
• I would advise against that, doing this accurately is nontrivial Jan 27, 2022 at 18:57
• Excellent point, I would crib the actual math from LAPACK itself (see ZLAEV2). (My point being, the OP might find it easier to crib that routine into their environment, than drill into / debug their whole stack) Jan 27, 2022 at 19:02