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()

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!?).

It looks like eig() instead chooses to make real, and positive, the vector entry with the largest abs value. 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 abs value avoids this.

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.