How to remove discontinuities from complex angle of NumPy eigenvector components?

Question

I am using NumPy's linalg.eig on square matrices. My square matrices are a function of a 2D domain, and I am looking at its eigenvectors' complex angles along a parameterized circle on this domain. As long as the path I am considering is smooth, I expect the complex angles of each eigenvector's components to be smooth. However, for some cases, this is not the case with Python (although it is with other programming languages). For the parameter M=0 (some argument in my matrix that appears on its diagonal), I have components that look like:

when they should ideally look like (M=0.1):

What I have tried:

• I verified that the matrices are Hermitian in both cases.
• When I use linalg.eigh, M=0.1 becomes discontinuous while M=0 sometimes becomes continuous.
• Using np.unwrap did nothing.
• The difference between component phases (i.e. np.angle(v1-v2) for eigenvector v=[[v1],[v2]]) is smooth/continuous, but this is not what I want.
• Fixing the NumPy seed before solving did nothing for different values of the seed. For example: np.random.seed(1).

What else can I do? I am trying to use Sympy's eigenvects just because I am running out of options, and I asked another question asking about another potential approach here: How do I force first component of NumPy eigenvectors to be real? . But, I do not know what else I can try.

Here is a minimal working example that works nicely in a Jupyter notebook:

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

M = 0.01; # nonzero M is okay
M = 0.0; # M=0 causes problems

def matrix_generator(kx,ky,M):
    a = 2.46;    t = 1;    k = np.array((kx,ky));
    d1 = (a/2)*np.array((1,np.sqrt(3)));d2 = (a/2)*np.array((1,-np.sqrt(3)));d3 = -a*np.array((1,0));
    sx = np.matrix([[0,1],[1,0]]);sy = np.matrix([[0,-1j],[1j,0]]);sz = np.matrix([[1,0],[0,-1]]);
    hx = np.cos(k@d1)+np.cos(k@d2)+np.cos(k@d3);hy = np.sin(k@d1)+np.sin(k@d2)+np.sin(k@d3);
    return -t*(hx*sx - hy*sy + M*sz)

n_segs = 200; #number of segments in (kx,ky) loop
evecs_along_loop = np.zeros((n_segs,2,2),dtype=float)
# parameterize circular loop
kx0 = 0.5; ky0 = 1; r1=0.2; r2=0.2; 
a = np.linspace(0.0, 2*np.pi, num=n_segs+2)
kloop=np.zeros((n_segs+2,2))
for i in range(n_segs+2):
    kloop[i,:]=np.array([kx0 + r1*np.cos(a[i]), ky0 + r2*np.sin(a[i])]) 
    
# assign eigenvector complex angles
for j in np.arange(n_segs):
    np.random.seed(2)
    H = matrix_generator(kloop[j][0],kloop[j][1],M)
    eval0, psi0 = LA.eig(H)
    evecs_along_loop[j,:,:] = np.angle(psi0)
    
# plot eigenvector complex angles
for p in np.arange(2):
    for q in np.arange(2):
        print(f"Phase for eigenvector element {p},{q}:")
        fig = plt.figure()
        ax = plt.axes()
        ax.plot((evecs_along_loop[:,p,q]))
        plt.show()

Clarification for anon01's comment: For M=0, a sample matrix at some value of (kx,ky) would look like:

a = np.matrix([[0.+0.j, 0.99286437+1.03026667j],
 [0.99286437-1.03026667j, 0.+0.j]])

For M =/= 0, the diagonal will be non-zero (but real).

Answer 1

I think that in general this is a tough problem. The fundamental issue is that eigenvectors (unlike eigenvalues) are not unambiguously defined. An eigenvector v of M with eigenvalue c is any non-zero vector for which

M*v = c*v

In particular for any non zero scalar s, multiplying an eigenvector by s yields an eigenvector, and even if you demand (as usual) that eigenvectors have length 1, we are still free to multiply by any scalar of absolute value 1. Even worse, if v1,..vd are orthogonal eigenvectors for c, then any non-zero linear combination of the v's is also an eigenvector for c.

Different eigendecomposition routines might well, therefore, come up with very different eigenvectors and still be doing their job. Moreover some routines might produce eigenvectors that are far apart for matrices that are close together.

A simple tractable case is where you know that all your eigenvalues are non-degenerate (i.e. each eigenspace is of dimension 1) and you happen to know that for a particular i, the i'th component of each eigenvector will be non zero. Then you could multiply the eigenvector v by a scalar, of absolute value 1, chosen so that after the multiplication v[i] is a positive real number. In C

s = conj(v[i])/cabs(v[i])
where 
conj(z) is the complex conjugate of the complex number z, 
and cabs(z) is the absolute value of the complex number z

Note that the above supposes that we are using the same index for every eigenvector, though the factor s varies from eigenvector to eigenvector.

This would impose a uniqueness on the eigenvectors, and, one would hope, mean that they varied continuously with the parameters of your matrix.