Eigenvectors of Lanczos Algorithm differ from numpy.linal.eig() for Complex Matrices

Question

I seek the lowest eigenvalues and eigenvectors for a complex-valued Hermitian matrix. For any real-valued Hermitian matrix, the following code works perfectly:

import numpy as np
import numpy.linalg as la
import copy as cp

#Generates the representation in Krylov subspace of a Hermitian NxN matrix using the Lanczos algorithm and an initial vector guess vg.
def Lanczos(H,vg):
    Lv=np.zeros((len(vg),len(vg)), dtype=complex) #Creates matrix for Lanczos vectors
    Hk=np.zeros((len(vg),len(vg)), dtype=complex) #Creates matrix for the Hamiltonian in Krylov subspace
    Lv[0]=vg/la.norm(vg) #Creates the first Lanczos vector as the normalized guess vector vg
     
    #Performs the first iteration step of the Lanczos algorithm
    w=np.dot(H,Lv[0]) 
    a=np.dot(np.conj(w),Lv[0])
    w=w-a*Lv[0]
    Hk[0,0]=a
     
    #Performs the iterative steps of the Lanczos algorithm
    for j in range(1,len(vg)):
        b=(np.dot(np.conj(w),np.transpose(w)))**0.5
        Lv[j]=w/b
         
        w=np.dot(H,Lv[j])
        a=np.dot(np.conj(w),Lv[j])
        w=w-a*Lv[j]-b*Lv[j-1]
        
        #Creates tridiagonal matrix Hk using a and b values
        Hk[j,j]=a
        Hk[j-1,j]=b
        Hk[j,j-1]=np.conj(b)
        
    return (Hk,Lv)

But for complex values, the output eigenvectors aren't the same (though the the eigenvalues are the same!):

H = np.random.rand(8,8) + np.random.rand(8,8)*1j #Generates random complex-valued 8x8 matrix
H = H + H.conj().T #Ensures 8x8 matrix is symmetric (hermitian)
Hc = cp.copy(H)
a,b = np.linalg.eig(H) #Directly computes eigenvalues of H
vg = np.random.rand(8) + np.random.rand(8)*1j #Generates random guess vector
Hk,Lv = Lanczos(Hc,vg) #Applies the Lanczos algorithm to H using the guess vector vg
A,B= np.linalg.eig(Hk)

#While the following outputs are the same
a[0]
A[0]
#The following are not!!!
b[:,0]
np.dot(B[:,0],Lv)

Any idea what I am doing wrong? Thank you.

--- Solution ---

Shoutout to Acccumulation for pointing it out, but it appears that there is nothing wrong in this procedure, as the eigenvalue check yields:

a[0]*b[:,0] - np.dot(H,b[:,0])
array([ 3.55271368e-15+1.11022302e-15j,  4.44089210e-16-3.65257395e-16j,
        0.00000000e+00+9.99200722e-16j, -4.44089210e-16+3.33066907e-16j,
       -2.22044605e-15-1.66533454e-16j,  1.33226763e-15+0.00000000e+00j,
        2.66453526e-15+1.22124533e-15j,  3.10862447e-15+1.22124533e-15j])
A[0]*np.dot(B[:,0],Lv) - np.dot(H,)
array([2.22044605e-15+1.77635684e-15j, 2.66453526e-15+1.55431223e-15j,
       2.66453526e-15+1.55431223e-15j, 2.66453526e-15+1.77635684e-15j,
       3.55271368e-15+2.66453526e-15j, 1.77635684e-15+1.99840144e-15j,
       1.33226763e-15+1.77635684e-15j, 8.88178420e-16+1.55431223e-15j])

implying that both b[:,0] and np.dot(B[:,0],Lv) are eigenvectors to the Hermitian matrix H.

Answer 1

It looks like this is the expected behavior. There is a function in numpy specifically for calculating the eigenvals/eigenvectors of complex Hermitian matrices: numpy.linalg.eigh(), so we wouldn't expect numpy.linalg.eig() to get the correct values for a Complex Hermitian matrix.

Answer 2

For a nondegenerate matrix, each eigenvalue has a one-dimensional space of eigenvectors associated with it. If the vector space has a norm, we can narrow our choice of what eigenvector to take by requiring that its norm be 1, but for vectors over the complex numbers, there are still an infinite number of eigenvectors with norm 1 for each eigenvalue. So it appears that the two algorithms choose different eigenvalues from among all the valid ones.

For degenerate matrices, the situation is even more complicated, as eigenvalues can have multidimensional eigenspaces.