Numpy.linalg.eig is giving different results than numpy.linalg.eigh for Hermitian matrices

Question

I have one hermitian matrix (specifically, a Hamiltonian). Though phase of a singe eigenvector can be arbitrary, the quantities I am calculating is physical (I reduced the code a bit keeping just the reproducible part). eig and eigh are giving very different results.

import numpy as np
import numpy.linalg as nlg
import matplotlib.pyplot as plt
   
def Ham(Ny, Nx, t, phi):
    h = np.zeros((Ny,Ny), dtype=complex)
    for ii in range(Ny-1):
        h[ii+1,ii] = t
    h[Ny-1,0] = t    
    h=h+np.transpose(np.conj(h))
    
    u = np.zeros((Ny,Ny), dtype=complex)
    for ii in range(Ny):
        u[ii,ii] = -t*np.exp(-2*np.pi*1j*phi*ii)
    u = u + 1e-10*np.eye(Ny)
    
    H = np.kron(np.eye(Nx,dtype=int),h) + np.kron(np.diag(np.ones(Nx-1), 1),u) + np.kron(np.diag(np.ones(Nx-1), -1),np.transpose(np.conj(u)))
    H[0:Ny,Ny*(Nx-1):Ny*Nx] = np.transpose(np.conj(u))    
    H[Ny*(Nx-1):Ny*Nx,0:Ny] = u                           

    x=[]; y=[]; 
    for jj in range (1,Nx+1):
        for ii in range (1,Ny+1):
            x.append(jj); y.append(ii)        
    x = np.asarray(x)
    y = np.asarray(y)
    return H, x, y

def C_num(Nx, Ny,  E, t, phi):
    H, x, y = Ham(Ny, Nx, t, phi)
    
    ifhermitian = np.allclose(H, np.transpose(np.conj(H)), rtol=1e-5, atol=1e-8)
    assert ifhermitian == True
    
    Hp = H
    V,wf = nlg.eigh(Hp)     ##Check. eig gives different result
    idx = np.argsort(np.real(V))
    wf = wf[:, idx]

    normmat = wf*np.conj(wf)
    norm = np.sqrt(np.sum(normmat, axis=0))
    wf = wf/(norm*np.sqrt(len(H)))
    
    wf = wf[:, V<=E]  ##Chose a subset of eigenvectors              
    
    V01 = wf*np.exp(1j*x)[:,None]; V12 = wf*np.exp(1j*y)[:,None]
    V23 = wf*np.exp(1j*x)[:,None]; V30 = wf*np.exp(1j*y)[:,None]
    
    wff = np.transpose(np.conj(wf))
    C01 = np.dot(wff,V01); C12 = np.dot(wff,V12); C23 = np.dot(wff,V23); C30 = np.dot(wff,V30)

    F = nlg.multi_dot([C01,C12,C23,C30])

    ifhermitian = np.allclose(F, np.transpose(np.conj(F)), rtol=1e-5, atol=1e-8)
    assert ifhermitian == True

    evals, efuns = nlg.eig(F)    ##Check eig gives different result
    C = (1/(2*np.pi))*np.sum(np.angle(evals));
    
    return  C

C = C_num(16, 16, 0, 1, 1/8)
print(C)
            













  

Changing both nlg.eigh to nlg.eig, or even changing only the last one, giving very different results.

Answer 1

As I mentioned elsewhere, the eigenvalue and eigenvector are not unique.

The only thing that is true is that for each eigenvalue $A v = lambda v$, the two matrices returned by eig and eigh describe those solutions, it is natural that eig inexact but approximate results.

You can see that both the solutions will triangularize your matrix in different ways

H, x, y = Ham(16, 16, 1, 1./8)
D, V = nlg.eig(H)
Dh, Vh = nlg.eigh(H)

Then

import matplotlib.pyplot as plt
plt.figure(figsize=(14, 7))
plt.subplot(121);
plt.imshow(abs(np.conj(Vh.T) @ H @ Vh))
plt.title('diagonalized with eigh')
plt.subplot(122);
plt.imshow(abs(np.conj(V.T) @ H @ V))
plt.title('diagonalized with eig')

Plots this

That both diagonalizations were successfull, but the eigenvalues are indifferent order. If you sort the eigenvalues you see they match

plt.plot(np.diag(np.real(np.conj(Vh.T) @ H @ Vh)))
plt.plot(np.diag(np.imag(np.conj(Vh.T) @ H @ Vh)))
plt.plot(np.sort(np.diag(np.real(np.conj(V.T) @ H @ V))))
plt.title('eigenvalues')
plt.legend(['real eigh', 'imag eigh', 'sorted real eig'], loc='upper left')

Since many eigenvalues are repeated, the eigenvector associated with a given eigenvalue is not unique as well, the only thing we can guarantee is that the eigenvectors for a given eigenvalue must span the same subspace.

The diagonalization test is the best in my opinion.

Is eigh always better than eig?

If you search for the eigenvalues in the lapack routines you will have many options. So it is I cannot discuss each possible implementation here. The common sense says that we can expect that the symmetric/hermitian routines to perform better, otherwise ther would be no reason to add one more routine that is more limited. But I never tested carefully the behavior of eig vs eigh.

To have an intuition compare the equation for tridiagonalization for symmetric matrices, and the equation for reduction of a general matrix to its Heisenberg form found here.