Eigenanalysis of complex hermitian matrix: different phase angles for EIG and EIGH

Question

I understand that eigenvectors are only defined up to a multiplicative constant. As far as I see all numpy algorithms (e.g. linalg.eig, linalg.eigh, linalg.svd) yield identical eigenvectors for real matrices, so apparently they use the same normalization. In the case of a complex matrix, however, the algorithms yield different results.

That is, the eigenvectors are the same up to a (complex) constant z. After some experimenting with eig and eigh I realised that eigh always sets the phase angle (defined as arctan(complex part/real part)) to 0 for the first component of each eigenvector whereas eig seems to start with some (arbitrary ?) non-zero phase angle.

Q: Is there a way to normalize the eigenvectors from eigh in the way eig is doing it (that is not to force phase angle = 0)?

Example

I have a complex hermitian matrix G for which I want to calculate the eigenvectors using the two following algorithms:

• numpy.linalg.eig for a real/complex square matrix
• numpy.linalg.eighfor a real symmetric/complex hermitian matrix (special case of 1.)

Check that G is hermitian

# check if a matrix is hermitian
def isHermitian(a, rtol=1e-05, atol=1e-08):
    return np.allclose(a, a.conjugate().T, rtol=rtol, atol=atol)

print('G is hermitian:', isHermitian(G))

Out:

G is hermitian: True

Perform eigenanalysis

# eigenvectors from EIG()
l1,u1 = np.linalg.eig(G)
idx = np.argsort(l1)[::-1]
l1,u1 = l1[idx].real,u1[:,idx]

# eigenvectors from EIGH()
l2,u2 = np.linalg.eigh(G)
idx = np.argsort(l2)[::-1]
l2,u2 = l2[idx],u2[:,idx]

Check eigenvalues

print('Eigenvalues')
print('eig\t:',l1[:3])
print('eigh\t:',l2[:3])

Out:

Eigenvalues
eig     : [2.55621629e+03 3.48520440e+00 3.16452447e-02]
eigh    : [2.55621629e+03 3.48520440e+00 3.16452447e-02]

Both methods yield the same eigenvectors.

Check eigenvectors

Now look at the eigenvectors (e.g. 3. eigenvector) , which differ by a constant factor z.

multFactors = u1[:,2]/u2[:,2]
if np.count_nonzero(multFactors[0] == multFactors):
    print("All multiplication factors are same:", multFactors[0])
else:
    print("Multiplication factors are different.")

Out:

All multiplication factors are same: (-0.8916113627685007+0.45280147727156245j)

Check phase angle

Now check the phase angle for the first component of the 3. eigenvector:

print('Phase angel (in PI) for first point:')
print('Eig\t:',np.arctan2(u1[0,2].imag,u1[0,2].real)/np.pi)
print('Eigh\t:',np.arctan2(u2[0,2].imag,u2[0,2].real)/np.pi)

Out:

Phase angel (in PI) for first point:
Eig     : 0.8504246311627189
Eigh    : 0.0

Code to reproduce figure

num = 2
fig = plt.figure()
gs = gridspec.GridSpec(2, 3) 
ax0 = plt.subplot(gs[0,0])
ax1 = plt.subplot(gs[1,0])
ax2 = plt.subplot(gs[0,1:])
ax3 = plt.subplot(gs[1,1:])
ax2r= ax2.twinx()
ax3r= ax3.twinx()
ax0.imshow(G.real,vmin=-30,vmax=30,cmap='RdGy')
ax1.imshow(G.imag,vmin=-30,vmax=30,cmap='RdGy')
ax2.plot(u1[:,num].real,label='eig')
ax2.plot((u2[:,num]).real,label='eigh')
ax3.plot(u1[:,num].imag,label='eig')
ax3.plot((u2[:,num]).imag,label='eigh')
for a in [ax0,ax1,ax2,ax3]:
    a.set_xticks([])
    a.set_yticks([])
ax0.set_title('Re(G)')
ax1.set_title('Im(G)')
ax2.set_title('Re('+str(num+1)+'. Eigenvector)')
ax3.set_title('Im('+str(num+1)+'. Eigenvector)')
ax2.legend(loc=0)
ax3.legend(loc=0)
fig.subplots_adjust(wspace=0, hspace=.2,top=.9) 
fig.suptitle('Eigenanalysis of Hermitian Matrix G',size=16)
plt.show()

Answer 1

As you say, the eigenvalue problem only fixes the eigenvectors up to a scalar x. Transforming an eigenvector v as v = v*x does not change its status as an eigenvector.

There is an "obvious" way to normalize the vectors (according to the euclidean inner product np.vdot(v1, v1)), but this only fixes the amplitude of the scalar, which can be complex.

Fixing the angle or "phase" is kind of arbitrary without further context. I tried out eigh() and indeed it just makes the first entry of the vector real (with an apparently random sign!?).

eig() instead chooses to make real the vector entry with the largest real part. For example, here is what I get for a random Hermitian matrix:

n = 10
H = 0.5*(X + X.conj().T)
np.max(la.eig(H)[1], axis=0)

# returns
array([0.57590624+0.j, 0.42672485+0.j, 0.51974879+0.j, 0.54500475+0.j,
       0.4644593 +0.j, 0.53492448+0.j, 0.44080532+0.j, 0.50544424+0.j,
       0.48589402+0.j, 0.43431733+0.j])

This is arguably more sensible, as just picking the first entry, like eigh() does, is not very robust if the first entry happens to be very small. Picking the max value avoids this. I am not sure if eig() also fixes the sign (a random matrix is not a very good test case for this as it would be very unusual for all entries in an eigenvector to have negative real parts, which is the only case in which an unfixed sign would show up).

In any case, I would not rely on the eigensolver using any particular way of fixing phases. It's not documented and so could, in principle, change in the future. Instead, fix the phases yourself, perhaps the same way eig() does it now.

Answer 2

In my experience (and there are many questions here to back this up), you NEVER want to use eig when eigh is an option - eig is very slow and very unstable. The relevance of this is that I believe your question is backward - you want to normalize the eigenvectors of eig to be like those of eigh, and this you know how to do.